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I .  INTRODUCTION 


The  BRLSC  code  is  a  modified  version  of  the  HELP 
code  and  was  developed  under  Contract  DAAD  05 -70-C-0275 ^ 
to  the  U.  S.  Army  Ballistics  Research  Laboratory  for  the  . 
purpose  of  predicting  phenomena  associated  with  the  perfor¬ 
mance  of  shaped  charges .  During  the  process  of  modifying 
HELP  to  enable  it  to  handle  shaped  charge  problems,  a  number 
of  additions  and  changes  were  introduced.  The  resulting 
differences  between  HELP  and  BRLSC  fall  into  two  basic 
categories.  The  f irs t . category  includes  the  necessary 
additions  to  enable  the  code  to  generate  shaped  charge 
problems,  treat  explosives,  simulate  a  slip  surface  between 
the  liner  and  the  explosive,  and  insure  realistic  jet  for-  . 
mation.  These  additions  were  specifically  required  by  the 
shaped  charge  application.  The  second  category  consists  of 
changes  made  to  improve  the. way  in  which  free  surface  cells 
are  handled,  to  allow  the  thin  liner  to  move  through  the 
'grid  without  becoming  dis torted , . to  increase  the  reliability 
and  stability  of  calculations,  and  to  decrease  .the  cost  of 
running  problems.  As  the  result  of  these  additions  and 
changes,  most  of  the  subroutines  in  BRLSC  have  been  changed 
substantially  from  the  original  HELP  version.  However,  since 
most  of  the  logic  involved  is  similar  to  that  employed  in 
HELP,  users' familiar  with  HELP  should  have  no  difficulties 
with  BRLSC. 
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It  is  the  purpose  of  the  present  document  (Volumes  I 
and  II)  to  report  the  BRLSC  program  in  detail  for  potential 
users.  A  companion  report  to  this  one,  by  R.T.  Sedgwick, 

J.M.  Walsh,  and  the  present  author,  gives  a  general  discussion 
of  the  shaped  charge  problem  and  computed  results  for  three 
specific  cases. 

In  Section  II  of  this  report  a  general  description  of 
the  numerical  method  is  presented.  Section  III  deals  with 
the  elastic-plastic  constitutive  relations.  In  Section  IV, 
the  method  used  to  propagate  tracer  particles  is  described 
and  the  way  in  which  these  particles  are  used  is  explained  in 
detail.  Section  V  describes  the  pressure  iteration  scheme 
used  in  determining  the  pressures  in  cells  containing  more 
than  one  material.  In.  Section  VI  the  method  used  to  apply 
a  pressure  boundary  condition  is  described.  (By  using  this 
pressure  boundary  condition,  an  effective  slip  surface  is 
introduced.)  Section  VII  describes,  several  sample  calculations 
which  have  been  run  with  BRLSC.  Finally,  Section  VIII  is  a 
users  guide  and  provides  a  brief  description  of  every  sub¬ 
routine,  detailed  instructions  for  generating  and  restarting 
problems,  and  a  dictionary  of  important  variables.  In  Volume 
II  of  this  report,  a  complete  listing  of  the  BRLSC  FORTRAN 
program  is  provided. 
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II.  GENERAL  DESCRIPTION  OF  THE  NUMERICAL  METHOD 


In  the  present  section  the  conservation  equations  and 
the  numerical  model  used  to  treat  these  equations  are  presented. 
In  addition  the  equations  of  state  used  to  model  both  the  con¬ 
densed  inert  material  and  the  explosives  are  presented,  along 
with  the  necessary -equation-of-state  parameters  for  ten  liner 
materials  and  eighteen  explosives. 


CONSERVATION  EQUATIONS 


Space  is  divided  into  fixed  Eulerian  cells  through 
which  the  fluid  moves.  To  arrive  at  expressions  for  the  rate 
of  change  of  total  mass,  momentum  and  energy  within  such  a 
cell,  it  is  convenient  to  start  with  the  equations  of  motion 
in  the  form: 


(a.  . u. ) 

\  ij  j) 


Here  the  are  the  components  of  the  stress  tensor,  which 

can  be  regarded  as  the  sum  of  the  hydrostatic  pressure,  -fi—P, 
and  the  deviatoric  stress  component,  s^j  ,  i.e., 


°ij  ■  sij  ■  5ijp 
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and  E.p  is  the  total  energy  (kinetic  plus  internal)  per  unit 
mass.  Tensor  notation  is  implied,  so  that  repeated  indices 
denote  summation. 

Expanding  the  convective  derivatives  in  Eqs .  (2)  and 
(3),  Df/Dt  =  3f/3t  +  u^  3f/3x^,  then  adding  Eq .  (1)  times 
to  Eq.  (2),  and  Eq.  (1)  times  E^  to  Eq.  (3),  and  collecting 
terms,  gives 


3 

/  \ 

3 

3  /  \ 

3 1 

(pUj) 

axi  aij 

axi  (puiuj) 

3 

3 1 

(pEp ) 

1  8x .  (puiEl) 

(5) 

(6) 


For  the  developments  to  follow  it  is  desirable  to  re¬ 
place  these  differential  equations  by  the  analogous  integral 
equations ,  obtained  by  integrating  over  the  cell  volume,  V, 
and  then  converting  the  volume  integral  of  divergences  to 
surface  integrals  over  the  cell  surfaces.  Equations  (1),  (5), 
and  (6)  then  become 


3 _ 

3t 


pdV 


puinids 


(7) 


3 

3 1 


3 

3 1 


pu.dV  = 
J 


pETdV  = 


s 


CTijnids  -  /  Puiujnid: 


aijujnids  -  f  puiETnids 


(8) 
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2,2  DIFFERENCE  EQUATIONS 


2.2.1  Division  into  Phases 

It  is  convenient  to  express  the  integral  conservation 
relations,  Eqs .  (7)  through  (9),  as  finite  difference  equations 
over  the  time  step  At  and  also  to  decompose  the  total  stress 
into  its  deviator  and  hydrostatic  components,  according 
to  Eq.  (4).  This  gives,  for  the  increments  of  total  mass  (m)  , 
momenta  ( mu ^ )  and  energy  (mE^)  within  the  cell 


s 


(12) 

Here,  the  terms  on  the  right  are  divided  into  increments  due 
to  pressure  and  stress  deviator  forces  on  the  cell  surface 
(first  colume)  ,  and  the  increments  (second  column)  due  to  the 
transports  of  mass,  momentum  or  energy  through  the  surface  of 
the  cell.  These  two  types  of  contributions  are  accounted  for 
in  jthe  computation  in  distinct  phases  .  Specifically,  during 
each  time  step  all  cells  are  updated  for: 

•  Pressure  and  deviator  stress  effects  in  Phase  1 

•  Transport  effects  in  Phase  2 
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Before  discussing  the  breakdown  of  the  calculations 
into  phases,  some  preliminary  definitions  will  be  useful. 
Superscript  n  on  a  variable  refers  to  the  value  of  the 
variable  at  the  beginning  of  the  time  step  and  superscript 
(n+1)  denotes  the  value  at  the  end.  In  this  discussion  a 
typical  cell  in  the  interior  of  the  grid  is  considered; 

Section  2. 2. 2. 6  discusses  the  special  conditions  which  de¬ 
scribe  the  handling  of  the  grid  boundaries  and  the  axis  of 
symmetry.  For  a  typical  cell,  denoted  by  a  value  of  the  index 
k,  the  dependent  variables  for  that  cell  are  written  P(k), 

Srr°°’  Szz^>  Srz(1°’  U(k)  >  v<k)  >  EmCk)  -  m(k)  >  P  00  , 
representing . respectively  the  hydrostatic  pressure,  the  three 

deviator  stress  terms,  the  radial  and  axial  components  of 
velocity,  the  specific  internal  energy,  the  mass,  and  the 
average  density  for  cell  k.  The  adjacent  cells  above,  below 
and  to  the  right  and  left  of  k  are  designated  respectively 
as  ka ,  kb,  kr ,  and  kl .  Here  the  terms  above,  below  right  and 
left  refer  to  a  cross-section  view  of  the  cells,  in  which  the 
left  border  is  parallel  to  the  axis  of  symmetry  with  z  in¬ 
creasing  upward  (see  Fig.  1).  Each  cell  is,  then,  the  torus 
obtained  by  rotating  the  rectangle  (since  Ar  f  Az  in  general) 
about  the  axis  of  symmetry.  In  the  discussion  that  follows, 
the  calculation  of  the  terms  on  the  right  of  Eqs .  (10)  through 
(12)  are  described  sequentially  starting  with  Phase  1. 

2.2.2  Phase  1- -Effects  of  Pressures  and  Deviator  Stresses 

In  Phase  1,  the  pressure  forces  and  stress  deviator 
forces  (if  effects  of  strength  are  being  included)  are  used,? 
to  update  the  velocities  and  internal  energies.  The  difference 
equations  used  in  Phase  1  were  chosen  because  they  prevent 
under-dense  cells  from  being  ’’kicked”  and  therefore,  eliminate 
the  need  for  "GLUEing".  An  artificial  viscosity  has  also  been 
included  in  the  differencing  scheme  in  order  to  minimize  insta¬ 
bility  in  low  velocity  stagnation  regions. 
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The  pressures  used  in  Phase  1  are  calculated  by  EQST 
and  set  in  subroutine  CDT.  The  deviator  stresses  used  are 
computed  in  subroutine  PH3.  If  the  effects  of  strength  are 
not  being  included,  the  deviator  stresses  are  set  to  zero. 

2.2.2. 1  Continuity  Equation,  (10) 

No  contribution  to  Phase  1. 


2.2. 2. 2  Equation  of  Motion,  (11) 

1.  Axial  Motion 

The  change  in  axial  momentum  is : 


A (mV)  = 


+  Qb 

-  sb  ) 

zz' 

-  a^p*  +  Q*  -  s\z) 

r 

At 

rz 

rz 

2.  Radial  Motion 

The  change  in  radial  momentum  is: 


A (mu)  = 


A 


(pA  + 

Q*  -  ) 

X  IT' 

-  Ar(pr  *  -  slt) 

AtSt 

-  AbSb  + 

(»r  -  A1) 

rz 

rz 

(POO 

-  srr(k)  + 

S,zW)]“ 
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2.2. 2.3  Energy  Equation 

The  change  in  total  cell  energy  is: 


“  + 

Q 

'  -  su  - 

ITS" 

zz  / 

rzj 

(pt 

+ 

Q*  -  S1  ) 

-  uV 

x  ZZ' 

rz 

( 

+ 

Q*  -  Sl  ) 

-  vlsl 

x  rr  t 

rz 

(pr 

+ 

Qr  -  S\. ) 

-  Vrs^ 

x  rr/ 

r  z 

2 . 2 . 2 . 4  Definition  of  Variables 

Am  -  area  of  cell  surface  m. 

P171  =  ’  pressure  at  the  interface  between  cell  k  and 
pm  =  [ p (k) P (km)  +  p(km)P(k)l/[p(k)  +  p(km)| 


sm  ,  sm 

rr 9  zz 


deviator  stress  at  the  interface  between  cell 
k  and  km. 

Szz  =  [p(k)S2Z(km)  +  p(km)Szz(k)]/[p(k)  +  p(km)j 

shear  stress  at  the  interface  between  cell  k 
and  km. 

S™z  =  |p(k)Srz(km)  +  p(km)Sr2(k)  /  p(k)  +  p (km) 
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vm,  um 


axial  and  radial  velocities  at  the  interface 
between  cell  k  and  km. 


V 


.m 


p (k) V (k)  +  p (km) V(km) 


/ 


p(k)  +  p (km) 


Qm  = 


pseudo-viscous  pressure  at  the  interface 
between  cell  k  and  km. 


Qm  =  C  *  [  U  (k)  -  U  (km) 
in.  radial  direction 

=  O  V(k)  -  V (km) 

in  axial  direction 


Jl 

(Pm-Sm  ) 

1 

'  rr/ 

p(k)  +  p(kn) 


/  2 


Pm-Sra  ) 
zzl 


p(k)  +  p  (km) 


/  2 


and  C  is  a  constant  =  1  if  interface  is  compressing 

=.4  if  interface  is  expanding 


2. 2. 2. 5  Extension  of  Phase  1  to  Mixed  Cells 


If  a  cell  contains  more  than  one  material,  it  is  necessary 
to  update  the  specific  internal  energies  of  each  material  in  the 
cell.  This  is  done  by  first  estimating  the  sound  speed  (c^) 
in  material  i.  Then 


ASIEi 


m(k)  *  AAIX 


m. 


pici 


where  AAIX  is  the  change  in  specific  internal  energy  for 
cell  k,  Is  the  mass  of  material  i,  and  is  the  density 

of  material  i. 
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The  densities  of  free  surface  mixed  cells  are  also  up¬ 
dated  in  Phase  1.  This  is  done  by  first  determining,  the  change 
in  volume  for  the  entire  cell  as  a  function  of  the  velocity 
field  according  to  the  formula 


AVOL  - 


aV 


-  ArUr 


AbVb 


-  aW1 


At 


The  density  of  material  i  is  then  updated  according  to 
the  following  expression: 


(n+1)  _ 


AVOL 


Zmi/pi 


n 


VCELL  p 


n2c?  rji 

1  1  JA-J  n2 


Pi  c? 


where  VCELL  =  total  cell  volume.  For  a  cell  containing  only 
one  material,  this  reduces  to 


(n+1)  _  Mi _ 

i  “  (1  +  AVOL/VCELL) 


2, 2. 2. 6  Special  Considerations  Along  the  Grid  Boundaries 

It  is  necessary  to  define  interface  values  for  stresses 
and  velocities  at  the  top,  bottom,  right  and  left  of  the  com¬ 
putational  grid. 

1.  Top  of  Grid  (transmi ttive) 

If  cell  k  is  at  the  top  edge  of  the  grid, 

P1  =  P(k)  Sb  =  S  _ (k) 

K  rz  rz'-  J 

sL  =  szzw  vt  =  v<k> 

UZ  =  U(k)  Q*  =  0 
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and  the  total  theoretical  energy  is  updated  by 


ETH  =  ETH  +  Ab  V1  (pt  -  sM  -  At 

L  '  zz f  rz 

2.  Right  of  Grid  (transmi ttive) 

If  cell  k  is  at  the  right  edge  of  the  grid 

Pr  =  P(k)  s'z  =  srz(k) 

S*  =  s  (k)  VT  =  V (k) 

rr  rrv  ^ 

UT  =  U(k)  Qr  =  0 

and  the  total  theoretical  energy  is  updated  by 

ETH  =  ETH  +  Ar|uT(pr  -  S*r)  -  VrS*z  At 

3.  Bottom  of  Grid  (transmi ttive) 

If  cell  k  is  at  the  bottom  edge  of  the  grid  and  the 
bottom  of  the  grid  is  transmittive  (CVIS  <  0) 

pb  =  P(kD  sbz  =  srz(k) 

Szz  =  SZZW  .  yb  =  V« 

Ub  =  U(k)  Qb  =  0 

and  the  total  theoretical  energy  is  updated  by 

ETH  =  ETH  -  Ab  Vb(pb  -  Sb  )  -  UbSb  At 

'  zz 1  rz 
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4.  Bottom  of  Grid  (reflective) 

If  cell  k  is  at  the  bottom  edge  of  the  grid  and  if 
the  bottom  of  the  grid  is  reflective  (CVIS  »  0) 


There  is  no  need  to  update  the  total  theoretical  energy  for  this 
case  since  no  work  is  done  across  a  reflective  grid  boundary. 


5.  Left  Edge  of  Grid  (Axis  of  Symmetry- -ref lective) 
If  cell  k  is  at  the  left  edge  of  the  grid 


There  is  no  need  to  update  the  total  theoretical  energy. 

2. 2.2.7  Additional  Modifications 

For  the  inadequately  resolved  region  in  the  immediate 
vicinity  of  the  axis,  a  change  was  made  to  prevent  unrealistic 
positive  radial  velocities  in  the  jet.  Specifically,  when  the 
jet  material  is  sufficiently  expanded  that  cell  pressures  are 
zero,  any  remaining  positive  radial  velocity  is  set  to  zero 
and  the  resulting  reduction  in  kinetic  energy  is  compensated 
for  by  increasing  the  internal  energy. 
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The  small  radial  velocities  which  are  corrected  by  this 
procedure  are  believed  to  be  caused  by  an  artificially  over¬ 
heated  jet,  the  artificial  heating  being  in  turn  traceable  to 
error  in  the  Phase  2  transport,  which  is  accurate  to  first 
order  only.  It  is  likely  that  the  correction  would  not  be 
needed  in  a  very  highly  resolved  calculation. 

2.2.3  Phase  2- -The  Effects  of  Transport 

In  Phase  2  the  transport  of  mass  ,  momentum  and  energy 
throughout  the  computational  grid  is  calculated.  This  trans¬ 
port  is  determined  by  integrating  the  last  terms  of  Eqs .  (10) 
through  (12).  In  the  discussion  below,  primary  attention  is 
given  to  the  case  of  transport  between  pure  cells.  At  the 
conclusion  of  this  discussion,  the  special  provisions  required 
to  treat  mixed  cells  are  described. 

2. 2. 3.1  Continuity  Equation 

The  transport  mass  is 

At 

s 

and  is  determined  for  each  cell  face  from 

6m  =  -  u^A^At  . 

is  the  density  of  the  cell  from  which  the  mass  moves  (donor 
cell)  ,  A^  is  the  area  of  the  face  and  is  an  interpolated 

value  of  the  velocity  component  normal  to  A^  representing 
approximately  the  velocity  at  the  interface  at  the  end  of  the 
time  step.  For  example,  considering  cells  k  and  ka 


/ 


puinjLds 


A  O)  =  - 
2 
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in  Fig.  1, 


1 

2 

vn(k)  +  vn(ka) 

1  + 

vn(ka). -  vn(k)  , J 

X  ' 

Az 

—  LA  L-j 

Calculated  transport  masses  are  subtracted  from  the  donor  cell 
mass  and  added  to  the  acceptor  cell  mass.  This  updating  is  • 
done  after  the  transport  .  terms  have  been  calculated;  all  trans 
port  terms  are  computed  using  the  pre - transport  quantities. 


2. 2.3. 2  Equation  of  Motion 
1.  Axial  Motion 

The  term  in  Eq.  (11)  for  axial  momentum  transport  is 

A  ( mu . 

2'  } 

5 

At  each  face  of  the  cell  the  transport  specific  momen¬ 
tum,  Uj ,  is  taken  to  be  the  axial  velocity  of  the  cell  from 
which  the  mass  moves  (donor  cell,  index  kd) ,  i.e.. 


-  6 1. 


pu .u- ) n . ds 
H  j  1/  l 


Uj  =  vn(kd)  „ 

Since  the  various  faces  of  the  cell  have  different  don6r  cells 
it  is  convenient  to  express  the  momentum  transport  for  each 
face 

6 (mv)  =  vn(kd)6m 


where 


is  the  mass  which  is  transported  across  the  cell  face,  as  given 
in  the  preceding  mass  transport  calculation.  Note  that  6 (mv) 
is  the  axial  momentum  and  the  6m  is  the  mass  transported 
in  either  the  r  or  z  direction,  depending  on  which  face  of 
the  cell  is  being  computed. 

2 .  Radial  Motion 

A2(mUj)  =  ‘  At  f  ( Puiu j)  ni(^s 

and,  by  analogy  with  the  axial  case,  the  equation  for  the 
transport  of  radial  motion  across  an  interface  is 

6 (mu)  =  un(kd)  6m 

where  un(kd)  is  the  time  n  velocity  of  the  donor  cell  and 
6m  =  -p^  u^A^At  is  the  mass  which  is  transported  across  the 
cell  face  in  question,  as  computed  in  the  continuity  equation 
above. 


2. 2. 3. 3  Energy  Equation 

From  Eq.  (12)  the  expression  for  transport  of  energy 

is 

A  (mE^j  =  -  At  f  (pu^E^jn^ds  . 

Vs 

To  evaluate  this  integral,  the  transported  specific  energy 
is  taken  to  be  that  of  the  donor  cell,  kd,  i.e., 


Ed  =  (kd) 


un (kd) 
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and  the  total  energy  transported  across  a  given  interface  is 
therefore  the  product  of  this  specific  energy  and  the  associated 
transported  mass  computed  from  the  continuity  equation,  i.e., 

A3(et)=  ED6m  . 

Once  the  amount  of  mass,  momentum  and  energy  to  be  transported 
is  known  for  all  faces  of  the  cell,,  the  cell  quantities  can  be 
updated  to  account*for  these'  effects .  New  cell  velocities  are 
then  determined  by  dividing  the  new  momenta  by  the  updated 
cell  mass.  This  fixes  the  new  kinetic  energy  so  that  the 
internal  energy  can  be  easily  determined  since  the  total  energy 
is  known. 

2.2. 3. 4  Extension  of  Phase  2  to  Mixed  Cells 

The  procedure  used  to  treat  mixed  cells  is  an  extension 
of  the  preceding  method  for  transport  between  pure  cells. 
Specifically,  the  mass,  say  for  material  b,  transported  across 
a  given  interface  is,  in  analogy  with  the  result  in  Section 
2. 2. 3.1  above,  just 

Pb  ui  ni  ds 
s 

where  new  is  the  donor  cell  density  for  material  b  and 

the  surface  integral  is  now  over  just  that  part  of  the  surface 
area  belonging  to  material  b.  The  density,  ,  is  determined 
by  a  pressure  iteration  which  is  the  subject  of  Section  V.  The 
determination  of  the  fractional  surface  area  belonging  to  a 
given  material  is  discussed  in  Section  IV,  paragraph  4.4. 

The  momentum  transport,  in  analogy  with  the  pure  cel! 
case,  is  simply  the  product  of  the  transported  mass  and  the 
appropriate  velocity  component  (axial  or  radial)  of  the  donor 
cell . 


A 


W  = At  f 
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The  energy  transported  is  the  product  of  the  mass 
transported  and  the  energy  per  unit  mass,  kinetic  plus  internal. 
The  kinetic-energy  per  gram  is  the  same  for  all  materials  in 
the  donor  cell  since  only  one  cell  velocity  is  computed,  but 
the  internal  energies  are  generally  different.  The  internal 
energies  are  updated  in  Phase  1  as  described  in  previous 
sections . 

As  in  the  case  of  pure  cells,  final  values  of  cell 
variables  in  the  mixed  cells  are  easily  determined  once  all 
the  transport  quantities  have  been  determined  for  the  cell. 

The  new  cell  velocity  components  are  calculated  by  dividing 
the  new  momenta  by  the  updated  cell  mass.  This  determines 
the  new  kinetic  energy  for  each  material.  The  internal 
energy  for  each  material  is  that  which  is  necessary  to  con¬ 
serve  total  energy,,  kinetic  plus  '  internal ,  for  the  individual 
materials. 


2 . 2 • 3 . 5  Additional  Considerations  for  Free  Surface  Cells 

If  a  mixed  cell  is  also  a  free  surface  cell,  an  addi¬ 
tional  flux  term  is  computed  for  each  of  the  interfaces  of  the 
cell.  This  is  an  associated  volume. flux  and  is  determined  by 
dividing  the  mass  flux  across  each  interface  by  the  donor  cell 
density  used  to  compute  the  nass  flux  across  the  interface. 

Thus,  for  each  interface  i,  the  volume  flux  for  material  n  is 

•  5V0L  =  -u.  A  At 
n.  l  n- 


where  A  is  the  fractional  area  of  material  n  at  inter- 
ni 

face  i.  P1I2  then  updates  the  density  of  material  n  in  a 
free  -surface  cell  by 


m  + 

,(n+l)  =  _ l 

n 


E  6m 


n  i. 


J!  +  E  6V0L 


n 

) 

n 


ni 
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2.3  EQUATIONS  OF  STATE 


There  are  two  equations  of  state  used  in  BRLSC.  The 
first,  used  for  determining  pressures  and  constant  energy 
compressibilities  in  the  liner,  is  that  due  to  J.  H.  Tillotson.^ 
Only  the  condensed  form  of  the  equation  of  state  is  included 
since  the  liner  material  should  not  be  shocked  enough  to  cause 
significant  vaporization  of  material.  The  second  equation  of 
state  is  for  the  explosive.  This  is  a  gamma-law  gas  equation 
of  state  modified  to  give  zero  pressures  ahead  of  the  detona¬ 
tion  front. 


2.3.1  Equation  of  State  for  the  Liner 


The  subroutine  EQST  contains  the  equation  of  state  for 

A 

the  liner  and  has  the  form,  for  r\  >  1  -  ; 


P  -  Pc  = 


a  + 


I  n; 
0 


+  1 


Ip  +  Ap  +  Bp2  , 


(13) 


and  for  the  constant-energy  compressibility 


(14) 


(15) 
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and 


In  these  equations  P,  I  and  p  are  pressure,  specific  inter¬ 
nal  energy  and  mass  density  respectively,  n  =  P/PQ  =  V  +  1, 

and  a,  b,  I  ,  A,  B,  and  p  are  constants  for  the  particular 
o  o 

material.  The  values  of  these  constants  for  ten  materials  are 
given  in  Table  1. 


It  is  necessary  to  use  the  preceedings  equations  for 
pressure  in  conjunction  with  those  for  compressibility  to 
insure  that  the  pressure  increases  monotonically  with  density, 
i.e : ,  (3P/3p) j  >  0.  This  is  necessary  to  insure  that  the  pressure 
iteration  scheme  for  calculating  pressures  in  mixed  cells 
converges . 

2.3.2  Equation  of  State  for  the  Explosive 


The  pressures  and  constant  energy  compressibilities  for 
the  explosive  are  computed  in  subroutine  EOSXPL  using  a  modified 
gamma-law  equation  of  state. 


P 


P 

ex 


(Y-I)P  •  MAX 


I 


min  ■ 


(i  - 


(17) 


and 


(18) 
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Dienes,  J.  K.,  et  al . ,  "An  Eulerian  Method  for  Calculating  Strength 
Dependent  Deformation,"  General  Atomic  Report  GAMD-8497,  February  1968. 


In  these  equations,  P,  I  and  p  are  the  pressure,  specific 
internal  energy,  and  mass  density  respectively,  and  y  and 
1^  are  constants  for  the  particular  explosive.  Values  of 
these  constants  for  eighteen  explosives  are  given  in  Table  2. 

*min*  an  -*-nPut  parameter  (usually  about  106  ergs) 
which  insures  that  (3P/3p)j  >  0,  The  variable  C  is  a  flag 
ranging  between  0  and  1  which  is  used  to  determine  whether  a 
cell  is  behind  (C  =  0) ,  intercepted  by  (0  <  C  <  1)  or  ahead 
(C  -  1)  of  the  detonation  front  (see  Fig.  2). 

Each  cycle,  the  position  of  the  detonation  front  is 
determined  by 


R ,  =  R  +  D  •  T 
d  o 


where  R,  is  the  radius  of  the  detonation  front,  R  is  the 
initial  radius  of  detonation,  D  is  the  detonation  velocity 
given  in  Table  2,  and  T  is  time.  For  each  cell  in  the  grid, 
Rmin>  t*ie  distance  fr0m  the  center  of  the  detonation  to  the 
corner  of  the  cell  k  closest  to  the  center  of  the  detona¬ 
tion  and  Rm„  .  the  distance  from  the  center  of  detonation 
to  the  corner  of  the  cell  k  furthest  from  the  center  of  detona¬ 
tion  are  calculated.  Then 


( 

r  /  R,  -  R  .  \ 

) 

<  1 ,  min 

0  d  min  ] 

U  ’  l  R  -  R  •  / 

( 

L  \  max  mm  /  _ 

) 

(19) 


By  using  this  technique  and  noting  that  cells  containing  ex¬ 
plosive  are  originally  set  up  with  1*1  and  p  =  p 

o  o 

the  pressure  in  a  cell  ahead  of  the  detonation  front  is  then 

P  =  (y-l)p  *  I  •  .  Since  there  is  a  pressure  cutoff  (Pmin) 

used  in  the  calculation  which  is  computed  to  be  2(y-l)p  *  Imin, 

o 

the  pressure  in  cells  ahead  of  the  detonation  are  set  to  zero. 
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TABLE  2 


EQUATION  OF  STATE  CONSTANTS  FOR  EXPLOSIVES1’ 


Material 

Code 

Number 

Explosive 

D 

cm/psec 

E 

0 

ergs/g 

P 

0 

g/ cm 

Y 

1 

Baratol 

.487 

1.108  x  10 10 

2.610 

3.42 

.2 

Comp  B,  Grade  A 

.798 

5.034 

1.717 

2.71 

3 

Comp  B-3 

.770 

5.423 

1.715 

2.54 

4 

Cyclotol  (77/23) 

.825 

5.271 

1.754 

2.73 

5 

Datb 

.752 

3.856 

1.780 

2.89 

6 

HMX 

.911 

6.395 

1.891 

2.74 

7 

LX-01-0 

.684 

3.087 

1.310 

2.93 

8 

LX- 04-1 

.847 

5.622 

1.865 

2.72 

9 

LX- 07-0 

.864 

5.627  ' 

1.865 

2.76 

10 

Nitroglycerine 

.770 

4.519 

1.600 

2.75 

11 

Nitromethane 

.629 

5.238 

1.128 

2.18 

12 

Octol  (77.6/22.4) 

.848 

5.134 

1.821 

2.83 

13 

PBX-9010 

.837 

5.283 

1.783 

2.76 

14 

PBX-9011 

.850 

5.453 

1.770 

2.76 

15 

PBX-9404-03 

.880 

6.409 

1.840 

2.65 

16 

PETIV 

.830 

4.993 

1.770 

2.81 

17 

RDX 

.864 

5.027 

1.767 

2.90 

18 

TNT 

.693 

3.729 

1.630 

2.73 

TLee,  E.  L.  ,  H.  C.  Hornig,  and  J.  W.  Kury,  "Adiabatic 
Expansion  of  High  Explosive  Detonation  Products," 
Lawrence  Radiation  Laboratory  Report  UCRL-50422,  May  2, 
1962. 


Center  of 
Detonation 


Position  of  Detonation 
at  Time  T 


Fig.  2 - -Calculation  of  C 


III.  NUMERICAL  METHOD  FOR  COMPUTING  DEVIATOR  STRESSES 

3 . 1  INTRODUCTORY  *  REMARKS 


If  it  is  desired  to  include  the  effects  of  material 
strength  in  a  particular  solution,  the  deviator  stresses 
used  in  Phase  1  in  conjunction  with  the  pressure  to  update 
the  velocities  and  internal  energies  must  be  recomputed  each 
cycle.  This  is  done  in  Phase  3  and  is  explained  in  the 
following  paragraphs. 
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3.2  DETERMINATION  OF  STRAIN  RATE  DEVIATORS 


As  a  first  step  in  computing  stress  deviators,  the 
instantaneous  strain-rate  deviators  are  computed  from  the 
velocity  field.  In  cylindrical  coordinates  (in  rectangular 
coordinates,  the  u/x  term  is  zero)  these  are 

e  =  v  -  i  ( u  +  v  +  —\ 
zz  y  3  \  x  y  x  / 

'rr  =  ux  *  I  (ux  +  vy  +  l) 

£rz  =  ( uy  +  vx)/2 

Here  the  cell  centered  velocity  gradients  are  computed  as  the 
difference  between  interface  velocities  (a  density  weighted 
average  velocity)  divided  by  the  cell  dimension.  For  example, 
the  interface  velocity  between  cells  k  and  kn  is  given  by 

n  =  P  (k)  v  (k)  +  p  (kn)  v  (kn) 
v  p00“+"p(k5j 
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The  above  calculation  provides  the  strain  rate  deviator 
associated  with  material  at  the  cell  center  at  the  beginning 
of  the  time  step.  It  is  necessary  to  recognize  that  .consti¬ 
tutive  equations  are  to  be  applied  only  along  a  particle 
path.  To  this  end,  one  must  take  a  convective  derivative. 
Consider  the  time  N  and  position  PN  of  a  particle  which, 
at  time  N+l,  will  be  at  the  cell  center.  At  time  N  this 
point  will  be  offset  from  the  cell  center  by 


6  x  =  -u(k)  At 
6y  =  -  v  (k)  At  . 


The  strain  rate  deviators  at  the  offset  point  are 

determined  by  interpolation  between  cell  center  strain  rate 
deviators.  This  is  done  by  weighting  the  contribution  of 
neighboring  cells  in  proportion  to  their  overlap  areas  with 
a  rectangle  of  cell  dimensions,  centered  at  PN.  For  example, 
in  the  above  sketch,  the  weighting  factors  are 

for  cell  k 


_ WK  =  (Ax  -  |6x|)(Ay  -  |«y|) 

for  cell  in  horizontal  direction 
WKH  =  | 5x | (Ay  -  | Sy | ) 
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for  cell  in  vertical  direction 
IV KV  =  ]  6y  ]  (Ax  -  |Sx|) 
for  cell  in  diagonal  direction 
WKD  =  | 6x | | 6y | 

so  that  the  resulting  interpolated  value  of  a  strain  rate 
deviator  component,  say  e  ,  is  given  by 

i. 


(WK)cr2(k)  +  (WKH)erz(kh)  +  (WKV)erz(kv)  +  (WKD)erz(kd) 

‘  WK"+"  WKH  +  IVKV_+  WKD 


Having  determined  ^zz>  err>  ^rz  f°r  ^ie  mass  which  will  be 
at  the  center  of  cell  k  at  the  end  of  the  time  step,  the 
deviatoric  strain  increments  can  now  be  calculated. 
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3.3  DETERMINATION  OF  STRESS  DEVIATORS 


The  deviatoric  strain  increments  are  next  used  in  the 
material  constitutive  equation,  to  update  the  stress  deviator 
tensor  sij •  T°  do  this  one  first  needs  to  compute  the  stress 
deviator  at  time  N  for  the  particle  in  question.  This  is 
done  by  employing  the  area  weighting  procedure  described 
above  for  strain  rate  deviators. 

The  deviatoric  stress  increments,  ds^.  ,  are  then  deter¬ 
mined  by  using  the  elastic  relation 


ds. -  =  2Gder. 
iJ 

where  G  is  the  modulus  of  rigidity  and  defj  are  the  incre¬ 
ments  of  deviatoric  strain.  When  such  an  increment  of  stress 
causes  the  yield  criterion  to  be  violated  each  stress  com¬ 
ponent  is  proportionately  reduced  to  bring  the  stress  state 
normally  back  to  the  yield  surface.  A  variable  yield  strength 

Y  =  (Ko  +  KiP)(l  -  E/Em) 


is  defined,  to  account  for  the  increase  in  strength  at  high 
pressures  and  the  decrease  in  strength  at  elevated  tempera¬ 
tures.  The  three  components  of  deviator  stress  are  updated 
by 


.N+l 

> 

zz 


N  , 
=  s  + 
zz 


2G  Ac 


zz 


N+l 

srr 

N+l 

s_ 


N 

s  +  2G  Ac 
rr  rr 


+  2GA£rz 
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The  following  yield  condition  is  imposed: 


s  .  .  s  .  -  =  s 2  +  s 2  +  s2 

1J  zz  rr  00 


+  2srz  i  2Y 


where  s00  =  -  (srr  +  szz)  anc*  ^  *s  t*ie  strength  in 

shear.  If  the  above  inequality  is  exceeded  during  a  time  step, 
as  would  occur  in  plastic  deformation,  all  of  the  deviatoric 
stress  components  are  proportionately  reduced.  Thus  the 
state  of  stress  has  been  returned  normally  to  the  yield 
surface.  If  the  second  invariant  is  smaller  than  2Y2,  the 
stress  state  is  in  the  elastic  regime  and  the  stresses  are 
not  reduced. 


The  preceding  calculation  updates  'the 
values  of  these  quantities  at  cell  centers  (as 
material  convection.)  These  updated  deviator 
then  used  along  with  the  hydrostatic  pressures 
recompute  the  velocities  and  internal  energies 


sij,  giving 
corrected  for 
stresses  are 
in  Phase  1  to 
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3.4  MODIFICATIONS  FOR  GRID  BOUNDARIES  AND  AXIS 


The  description  given  above  applies  to  normal  cells 
within  the  grid.  Additional  statements  must  be  made  to  cover 
grid  boundaries  and  the  axis  of  symmetry. 

Grid  boundaries  may  be  reflective  (bottom  grid  boundary) 
or  transmittive  (bottom,  right  or  top).  In  calculating  the 
strain  rate  deviators  at  reflective  grid  boundaries  (from  the 
velocity  field)  it  is  assumed  that  a  row  of  fictitious  cells 
exists  outside  the  grid.  Each  cell  in  this  row  is  assigned 
the  same  tangential  component  of  velocity  and  an  equal  and 
opposite  normal  component  of  velocity  as  its  real  neighbor 
cell  at  the  grid  boundary.  Thus,  for  cell  k  at  a  reflective 
bottom  grid  boundary,  the  fictitious  cell  just  below  k  has 
velocities 

u  =  u  (k) 
v  =  -v (k) 

These  velocities  are  then  used  in  an  otherwise  standard  (above) 
calculation  of  the  strain  rate  deviator  tensor  for  cell  k. 

For  transmittive  boundary  cells,  one  assumes  that  the. fic¬ 
titious  cell  has  both  of  its  velocity  components  equal  to 
the  velocity  of  the  cell  just  inside  the  boundary. 

u  -  u(k) 

v  =  v(k) 

With  these  assumptions,  it  is  then  possible  to  compute  strain 
rate  deviators  for  the  border  cells. 

The  axis  of  symmetry  is  essentially  a  reflective  grid 
boundary.  In  calculating  strain  rate  deviators  for  the 
column  of  cells  next  to  the  axis,  the  cell  ’'to  the  left”  is 
assumed  to  have  the  same  axial  velocity  and  an  equal  and 
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opposite  radial  velocity,  i.e., 


v  =  v(k) 
u  =  -u(k)  . 
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IV.  THE  PROPAGATION  OF  TRACER  PARTICLES 
AND  THEIR  USE  TO  DESCRIBE  MATERIAL  INTERFACES 


4.1  INTRODUCTORY  REMARKS 

The  most  complicated  logic  of  the  BRLSC  code  is  found 
* 

in  the  subroutines  dealing  with  the  definition  and  use  of 
material  interfaces.  Conceptually  the  task  is  not  difficult. 

The  idea  of  moving  tracer  particles  with  a  locally  defined 
velocity  and  using  them  to  describe  the  interface  between 
materials  is  a  simple  one.  Likewise,  employing  the  partial 
cell  boundary  areas  determined  by  interface  position  to  com¬ 
pute  the  mass  flux  of  materials  is  a  very  natural  use  of  these 
interfaces.  The  complexity  of  the  task  comes  with  making 
the  treatment  completely  general.  As  the  code  is  written, 
there  are  no  restrictions  on  the  direction  of  the  flow,  on 
the  shape  of  the  material  interface,  on  the  number  of  times 
an  interface  can  cross  a  cell  boundary,  or  on  the  number  of 
materials  contained  in  a  cell. 

It  is  hoped  that  the  present  section  will  aid  the  user 
in  understanding  the  organization  and  logic  of  these  routines. 

In  particular,  this  section  covers  the  definition  of  the  material 
interfaces,  the  computation  of  fractional  cell  boundary  areas 
for  each  material  in  a  mixed  cell,  and  the  use  of  these  frac¬ 
tional  areas  in  defining  the  mass  transport  terms.  The 
creation  of  mixed  and  pure  cells,  and  the  propagation  of  the 
tracer  particles  that  define  the  package^  boundaries  will  also 
be  discussed. 


INFACE,  NEWMIX,  NEKRHO,  FLUX,  CALFRC ,  ADJFLX ,  MOVTCR 

^"To  distinguish  the  INFACE  and  PH 2  (Phase  2)  functions,  it 
should  be  noted  that  the  mass  transports  (for  mixed  cells) 
are  only  computed  and  stored  in  INFACE.  The  actual  trans¬ 
port  of  all  quantities  ,  and  the  associated  updating  of  cell 
variables,  is  performed  in  PI  1 2  for  both  mixed  and  pure  cells. 
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4.2  DEFINING  MATERIAL  INTERFACES 


A  set  of  massless  tracer  particles  are  initially  posi¬ 
tioned  along  the  boundary  of  each  material  package.  These 
tracers  are  numbered  such  that  the  package  is  on  the  left  as 
one  proceeds  between  any  two  consecutive  tracers .  The 
material  interface  therefore  is  defined  by  the  line  segments 
connecting  these  tracers.  The  initial  density  of  tracers  is 
controlled  by  input  variables. 

To  prevent  the  boundaries  of  adjacent  packages  from 
crossing,  the  tracers  along  the  common  boundary  coincide 
exactly.  Because  these  identical  points  are  moved  with 
exactly  the  same  velocity,  they  remain  superimposed  during 
the  entire  calculation. 

Subroutine  INFACE  is  called  each  cycle  to  handle  these 
massless  tracer  particles  and  to  compute  the  effects  that 
the  position  of  these  particles  have  of  mass  fluxes,  etc. 
INFACE  itself  calls  a  number  of  subroutines  among  which  are: 
MOVTCR  which  moves  the  tracers  with  the  local  velocity  field, 
CALFRC  which  computes  the  fractional  cell  areas  defined  by 
the  position  of  the  tracers,  FLUX  which  computes  the  mass 
flux  across  mixed  cell  interfaces  and  ADJFLX  which  computes 
the  mass  fluxes  necessary  to  transport  all  the  mass  out  of  a 
cell  which  has  lost  an  interface. 
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4.3  ADVANCING  TRACER  PARTICLES 

The  first  thing  that  INFACE  does  each  subcycle  is 
to  call  MOVTCR  to  move  the  tracer  particles. 

Each  material ‘is  circumscribed  by  a  series  of  massless 
tracer  particles,  which  are  propagated  each  time  step  a 
distance1* 

Ax^  =  u.*At 

where  tu  is  a  local  average  velocity  vector  for  the  con¬ 
tinuum,  determined  by  an  area  overlap  method  which  gives  weight 
to  velocities  in  the  surrounding  cells.  Specifically,  a  rec¬ 
tangle  of  cell  dimensions  is  superimposed  on  the  particle  to 
be  moved  and  then 

E  w(L)  A  p(L) 

+  L  L 

ui  =  - ^ - 

E  A  p(L) 

L  L 


L  -  1,  . 4 


noise  . 
this 


Subroutine  IXFACE  is  subcycled  to  minimize  transport 

•4*  . 

Since  the  tracer  coordinates  are  in  grid  line  units, 
distance  is  converted  from  centimeters  to  cell  units 
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where  p(L)  is  the  density  and  w(L)  is  the  ce 11  -  centered 
* 

velocity  of  the  overlapped  cell  L  and  is  the  over¬ 

lap  area.  This  procedure  gives  a  spatially  continuous  velo¬ 
city  field  for  particle  propagation. 

The  passive,  cell-centered  tracers  (XP,  YP) ,  used  only 
for  plotting  the  mass  flow  within  the  package  boundaries ,  are 
moved  by  the  same  method. 

4.3.1  Additional  Modification  to  Insure  Realistic 

Jet  Definition  ~ 

An  additional  consideration  must  be  made,  when  moving 
the  interface  tracers,  to  insure  that  the  jet  is  well  defined. 
Since  the  tracers  are  moved  with  an  area  weighted  velocity 
and  since  a  very  large  velocity  gradient  exists  near  the 
stagnation  point  at  early  times,  some  of  the  tracers  describing 
the  inside  surface  of  the  liner  near  the  axis  aren't  moved 
downward  fast  enough.  This  error  arises  when  the  o\rerlap  area 
extends  into  the  relatively  low  velocity  region  near  the  stag¬ 
nation  point  and  causes  a  very  uneven  interface  to  form.  An 
example  of  this  can  be  seen  in  Fig.  3.  To  overcome  this 
problem,  a  special  modification  has  been  included  in  the  tracer 
moving  subroutine  (MOVTCR)  of  BRLSC.  This  is  outlined  below: 

1.  When  the  problem  is  first  set  up,  a  function 
describing  the  position  of  the  tracers  along 
the  inside  surface  of  the  liner  is  derived, 

(y  =  f (x) ,  x,  y  in  cm)  and  the  free  surface 
tracers  which  are  defined  by  this  function 
identified.  i.e.,  the  tracers  TX (NVOID, NS)  , 

TY  (NVOID ,  NS)  through  TX  (NVOID  ,.\D)  ,  TY  (NVOID  ,ND) 
define  the  inside  surface  of  the  tracer. 


These  ce 1 1 - centered  velocities  hare  been  updated  by  PHI 
(pressure  and  strength  effects)  for  this  time  step. 
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2.  In  MOVTCR,  the  physical  coordinates  (cm)  for 
each  tracer  are  computed.  Thus,  when  moving 
the  free  surface  tracers,  ,  y^  is  the  physi¬ 
cal  location  of  tracer  TX(NVOID,N),  TY(NVOID,N) 


3. 

4. 


The  tracers  are  then  moved  so  that  tracer  N 
moves  from  (xjj,  y*  )  to  (xj]+1,  yjj+1)  . 

yJJ+^  for  free  surface  tracers  in  the  range 
ND  N  >  NS  (defined  in  Step  1)  are  then 
repositioned 


=  MIN 


yn+1 


yn+1 

yN-l 


5.  The  new  physical  coordinates  for  the  tracers 
are  then  converted  -back  to  cell  coordinates. 


6.  The  coordinates  of  the  tracers  for  material 
package  1  (the  liner)  are  set  equal  to  the 
appropriate  free  surface  tracers  moved  with 
the  above  method. 
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4.4  COMPUTING  FRACTIONAL  AREAS 


4.4.1  •  Problem  Definition 

When  an  Eulerian  cell  is  cut  by  the  material  N  inter¬ 
face,  we  are  interested  in  knowing  what  fraction  of  each  side 
of  the  cell  should  be  associated  with  material  N. 


Eor  example,  the  fractional  areas  associated  with  material 
N  in  the  diagram  above  are  as  follows:  • 

Left  =  0 

Top  %  .  8tt  [x?  -  xi-i) 

Right  *  .7(yj  -  Yj-i)  2lT  xi 
Bottom  ^  .3tt^x?  -  x? 

These  fractional  areas  are  used  in  computing  the  mass  flux 
of  material  N  across  each  side  of  the  mixed  cell.  (Actually, 
for  each  mixed  cell  it  is  necessary  to  compute  and  store  only 
the  fractional  areas  for  the  right  and  top  sides  since  the 
cells  below  and  on  the  left  provide  the  information  for  the 
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other  two  sides.)  The  method  given  below  for  computing  these 
fractional  areas  puts  no  restrictions  on  how  many  interfaces 
are  in  a  cell  or  on  how  many  times  a  single  interface  crosses 
a  cell  boundary.  It  does,  however,  assume  that  a  material 
package  boundary  never  crosses  itself,  and  experience  has 
shown  that  when  it  does  the  logic  of  INFACE  breaks  down. 

This  is  not  a  limitation  of  the  method  since  the  boundary  will 
not  cross  itself  while  running  a  shaped  charge  problem. 

4.4.2  Computing  Intersections  of  Interfaces 

with  Grid  Lines 

During  each  subcycle  through  INFACE,  subroutine  CALFRC 
is  called  to  compute  fractional  cell  areas  for  mixed  cells. 

Subroutine  CALFRC  considers  consecutively  (two  at  a 
time)  the  tracer  points  which  circumscribe  each  material 
package  in  a  counter-clockwise  direction.  For  simplicity, 
the  coordinates  of  the  tracers  are  in  grid  line  units  rather 
than  centimeter  units.  Thus  the  coordinates  of  a  pair  of 
tracers  immediately  indicate  whether  they  straddle  one  or 
more  cell  boundaries. 

TX  and  TY  are  the  FORTRAN  variable  names  of  a  tracer’s 
x  and  y  coordinates,  respectively.  The  arrays  are  doubly 
dimensioned;  the  first  dimension  identifies  the  material 
package,  the  second  gives  the  sequential  ordering  of  the 
points.  In  the  example  shown  below.,  the  indices  indicate  that 
both  tracer  points  belong  to  material  package  one  (first  index) 
and  that  they  are  consecutive  points  (second  index),  namely 
the  20th  and  21st  points  describing  package  one.  The  x- 
coordinates  indicate  the  line  through  the  points  does  not 
intersect  a  vertical  grid  line,  whereas  the  y- coordinates 
indicate  the  line  through  the  points  intersects  the  fourth 
horizontal  grid  line. 
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TX ( 1 , 20)  =  2.5' 
TY(1 ,20)  =  3.5  . 

TX (1 , 21)  =  2.8] 
TY ( 1 , 2 1)  =  4.6  j 


4.4.3  Defining  Fractional  Areas  of  Intersected 

Cell  Boundaries 

Subroutine  CALFRC  defines  the  point  of  intersection 
of  the  line  between  two  tracers  and  a  grid  line  and  uses  it 
to  define  the  fractional  area  of  the  particular  cell  boundary 
associated  with  material  package  one.  The  code  recognizes 
from  the  sequential  numbering  of  the  points  that  material  one 
is  on  the  left  of  the  point  of  intersection. 

The  FORTRAN  variables  FRACTP  and  FRACRT  store  the 
fractional  areas  of  the  top  and  right  boundaries,  respectively, 
for  all  the  materials  of  a  mixed  cell.  These  arrays  are  also 
doubly  dimensioned  and  are  written  as  FRACTP (N,M)  and 
FRACRT (N,M)  where  the  first  dimension,  N ,  identifies  the 
material  package  number  and  the  second  dimension,  M,  is  an 
index  that  links  the  mixed  cell  K  to  the  special  mixed  cell 
arrays  by  the  relation,  M=FLAG (K) - 100 .  (See  Appendix  C.) 
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If  material  package  two  is  adjacent  to  material  package 
one,  the  material  two  tracers  are  identical  to  the  material 
one  tracers,  but  are  numbered  in  opposite  order.  For  example, 
package  two  tracers  (say  30  and  31)  that  coincide  with  the 
package  one  tracers  discussed  previously  would  be  as  follows: 

TX (2 , 30 )  =  2.8 
TY (2,30)  -  4.6 
TY (2 , 31)  =  2.5 
TY (2 , 31)  =  3.5 

The  fractional  area  to  the  right  of  the  intersection  point 
is  associated  with  material  package  two,  and  is  stored  in 
FRACRT  (2 ,M) . 

4.4.4  Defining  Fractional  Areas  of  Boundaries 
Not~  Intersected- 

Because  the  fractional  area  is  a  term  in  the  flux 
equations  applied  to  mixed  cells,  it  must  be  defined  for  the 
right  and  top  boundary  of  every  mixed  cell  even  if  one. or  both 
are  not  intersected  by  a  material  interface.  In  that  case  the 
fractional  area  either  equals  the  total  cell  boundary  area  or 

•  "(xi  ■  xi-i) 

=  h  ■  yj-i)2’  xi 

=  0. 

=  0. 


is  zero,  e.g. 


Material  Package  1 


Material 
Package  2 


FRACTP (1 ,M) 
FRACRT (1,M) 

FRACTP ( 2, M) 
FRACRT ( 2, M) 


42 


1&3  L_T^O  L_J®LJ  L_© 


L- 


n 

9 

y,,j 


“1 


§ 

~1 


In  the  above  example,  none  of  the  material  of  package  two  is 
transported  across  the  right  or  top  boundaries  of  the  cell. 

On  the  other  hand,  the  full  cell  area  is  used  to  compute  the 
mass  flux  of  the  material  of  package  one  across  these  boun¬ 
daries.  There  is  clearly  a  need  to  systematically  define 
the  fractional  areas  of  mixed  cell  boundaries  that  are  not 
intersected  by  a  material  interface.  The  following  discussion 
illustrates  how  the  fractional  area  of  these  non- intersected 
boundaries  are  correctly  defined  by  "presetting"  them  when 
the  interface  enters  a  cell  and  "resetting"  them  when  it 
leaves . 

The  direction  of  an  interface,  given  by  the  sequential 
order  of  the  tracer  points  ,  determines  whether  one  is  entering 
a  cell  or  leaving  it.  When  entering  a  cell,  and  if  crossing 
this  cell  boundary  for  the  first  time,  the  program  processes 
the  other  sides  of  the  cell  in  a  clockwise  order  and  presents 
their  fractional  areas  to  equal  the  total  cell  boundary  area 
until  it  comes  to  a  side  that  has  already  been  crossed  by  the 
same  material  interface  (i.e.,  a  side  which  has  an  associated 
fractional  area  that  is  non-zero  and  is  less  than  the  total 
area.) 

Fractional  areas  preset: 

Bottom  =  7t  (x?  -  xi -i 
Left  =  (y.  -  yj.1)21r 
Top  =  ir(x?  -  x?^) 


Case  (1) 
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Case  (2)  Fractional  areas  preset: 


If  the  interface  of  material  N  crosses  a  cell 
boundary  for  the  second  time,  the  code  does  not  preset  the 
fractional  areas  when  material  N  is  between  the  two  points 
of  intersection  (Case  4).  If  material  N  is  outside  the  two 
points,  the  fractional  areas  are  preset  as  described  above 
(Case  S) .  The  program  senses  that  material  N  is  between 
the  points  when  the  sum  of  the  fractional  area  computed  on 
the  previous  crossing (s)  and  the  one  computed  currently  is 
greater  than  the  total  area. 
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Actually  this  rule  holds  regardless  of  how  many  times  the 
interface  crosses  the  boundary,  as  illustrated  by  Case  6 
and  Case  7. 

Case  (6)  Fractional  areas  preset: 

Bottom  -  71  (xf  ‘  xi-i) 

.  Left  =  (y.  -  yj.1)2ir  x. 

Sum  of  areas 
<  Total  area 


Case  (7) 


Fractional  areas  preset: 
None 


When  the  interface  enters  cell  K  from  the  right, 
it  has  just  left  cell  K+l,  and  the  program  must  "reset"  to 
zero  the  fractional  areas  of  certain  non- intersected  sides 
of  ..cell  K+l , 
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In  Cases  8  through  10  the  side  where  the  interface 
leaves  is  crossed  only  once  by  this  interface.  Proceeding 
in  a  clockwise  direction  from  the  side  the  interface  has 
left,  the  program  resets  to  zero  the  fractional  areas  of  the 
non^intersected  sides  until  it  encounters  a  side  that  is 
intersected  by  this  interface. 


Case  (8) 


Fractional  areas  reset: 

T  op  =  0  . 

Right  -0. 


Fractional  areas  reset: 


Top  =  0. 


Fractional  areas  reset: 


Case  (IQ) 


In  Cases  11  and  12  below,  the  side  where  the  interface 
leaves  has  been  previously  crossed  by  this  interface.  The 
program  sums  the  fractional  areas  resulting  from  a  previous 
crossing  of  that  boundary  and  the  current  crossing.  If  the 
sum  of  these  areas  is  less  than  the  total  cell  boundary  area 
(Case  11)  none  of  the  fractional  areas  of  the  other  sides  are 
reset.  If  the  sum  is  greater  (Case  12),  the  program  proceds 
to  reset  to  zero  the  areas  of  the  other  non- intersected  sides 
of  the  cell. 

Case  (11)  Fractional  areas  reset:  . 

None 

revious 

Sum  of  areas 
<  Total  area 

revious 


Case  C12) 


Fractional  areas  reset: 
Top  =  0. 

Right  =  0  . 


The  calculation  of  the  fractional  areas  is  thus  re 
duced  to  considering  four  directions  of  an  interface  with 
respect  to  a  specified  cell  K: 

1.  entering  cell  K,  leaving  cell  KR 

2.  entering  cell  KR,  leaving  cell  K 

3.  entering  cell  K,  leaving  cell  KA 

4.  entering  cell  KA,  leaving  cell  K. 
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4.5  COMPUTING  THE  FLUX  TERMS  FOR  INTERFACE  CELLS 


The  computation  of  fractional  cell  boundary  areas  is 
performed  in  order  to  use  the  material  interface  positions 
to  define  the  mass  flux  of  materials  across  interface  cell 
boundaries.  For  pure  cells  the  mass  flux  (computed  in 
Phase  2)  between  two  cells  is 

Am  =  p,p  V,p  A  At 

where  is  the  donor  cell  density,  Vj  is  an  average  of 

the  cell-centered  velocities  of  the  two  cells,  and  A  is 
the  area  of  the  cell  boundary  for  which  the  mass  flux  is 
being  computed. 

For  cells  crossed  by  an  interface  this  equation 

becomes 

AmN  =  PTN  VT  FAN  At/S 

where  p^,^  is  the  donor  material  N  density,  FA^  is  a  frac 
tional  cell  boundary  area  for  material  N,  and  the  time  step 
At,  is  divided  by  S,  the  number  times  the  INFACE  subroutine 
is  subcycled.  Thus  an  interface  cell  has  mass  flux  terms 
(which  are  calculated  in  subroutine  FLUX)  for  its  right  and 
top  boundaries  associated  with  each  material  package  in  the 
grid. 
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Material  Package  2 


Material  Package  1 


f 


Void 


In  the  case  shown  above,  assume  there  are  two  material 
packages  in  the  grid.  There  are  four  mass  flux  terms  asso¬ 
ciated  with  cell  K.  Using  the  FORTRAN  variables  and  assuming 
M  is  the  location  in  the  mixed  cell  arrays  associated  with 
cell  K,  the  right  boundary  fluxes  would  be 

|  SAMMP (1  ,M)  |  ^  pT1  UT  .292  (y..  -  Yj_1)2.i\  At/s 

|  SAMMP  (2  ,M)  |  -v  pT2  UT  .  3 7 s(  y j  2tt  x±  At/s 

and  the  top  boundary  fluxes  would  be 

|  SAMPY  (1  ,M)  |  -v  pT1  VT  .  417  (x?  -  x?_1)ir  At/s 

|  SAiMPY (2  ,M)  |  a,  pT2  VT  .  S85 (x?  -  x?^)*  At/s 

The  sign  of  these  fluxes  indicates  the  direction  of  the  flow. 
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4.6  CREATING  MIXED  CELLS 


4- 


4.6.1  Assigning  Storate  Location 

Initially  those  cells  containing  a  material  interface 
are  created  in  SETUP  when  the  problem  is  generated.  Later, 
as  the  interfaces  move  across  the  grid,  new  mixed  cells  are 
created  in  NEWMIX  (called  by  CALFRC  which  is  called  by  INFACE) . 
Whenever  a  cell  boundary  is  straddled  by  a  pair  of  consecu¬ 
tive  tracers,  INFACE  checks  that  both  cells  (K  and  K^)  on 
either  side  of  the  boundary  are  mixed.  If  one  or  both  are 

pure  (i.e.,  MFLAG(K)  <  100,  or  MFLAG(K')  <  100),  subroutine 

* 

NEWMIX  is  called  and  an  unused  storage  location  in  the  mixed 
cell  arrays  is  assigned  to  that  cell.  [The  maximum  number  of 
mixed  cells  at  any  one  time  is  controlled  by  the  input  variable 
NMXCLS  which  should  correspond  to  the  dimensions  of  the  mixed 
cell  arrays.  If  the  user  tries  to  generate  more  than  NMXCLS 
mixed  cells,  NEWMIX  calls  for  an  error  exit.)  Before  a  cell 
K  becomes  mixed,  its  flag,  MFLAG(K),  indicates  what  material 
package  it  belongs  to.  NEWMIX  uses  that  information,  trans¬ 
ferred  through  a  variable  MO  in  blank  common,  to  define 
the  mixed  cell  variables  for  the  neitfly  mixed  cell,  as  illustrated 
by  the  following  FORTRAN  statements. 

MO  =  MFLAG(K)  (when  pure) 


(M  =  storage  location  for  cell  K  in  mixed  cell  arrays) 


MFLAG(K)  =  M+100 
XMAS S (MO, M)  =  AMX(K) 

SIE (MO ,M)  =  AIX(K) 

RHO (MO ,M)  =  AMX (K) / [TAU (I)*DY (J) ] 


(flag) 

(mass ) 

(specific  internal  energy) 
(density) 


A  storage  location  M  is  not  in  use  if  RIIO(l,M)  =  -1. 
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A  cell  that  becomes  "mixed"  in  INFACE  has  only  one  material 
until  subroutine  PH2  is  executed  and  the  other  materials  are 
transported  into  it.  Note  also  that  PH2  will  not  transport 
material  of  another  package  into  a  pure  cell.  The  interface 
of  the  other  package  must  first  cross  the  cell's  boundary  in 
INFACE  before  the  cell  can  receive  the  material  of  the  other 
package  in  PH2. 

4.6.2  Defining  Free  Surface  Cells 

A  cell  containing  a  free,  surface  often  contains  only 
one  material  and  yet  is  still  treated  as  a  mixed  cell  since 
an  interface  crosses  its  boundaries  and  the  fluxes  across  its 
boundaries  are  a  function  of  the  fractional  areas  computed  in 
INFACE  by  CALFRC.  Of  course,  it  is  possible  for  a  cell  con¬ 
taining  a.  void  (i.e.,  a  free  surface)  to  contain  more  than  one 
material.  A  "free  surface"  cell  K  is  flagged  by  setting 
RHO(NVOID,M)  =  1.,  where  NVOID  is  one  more  than  the  number 
of  material  packages  in  the  problem  and  M  =  MFLAG (K) - 100 . 

4.6.3  Accounting  for  Subcycles 

INFACE  usually  is  subcycled,  i.e.,  the  routine  is 
executed  several  times  each  cycle  using  a  fraction  of  the  time 
step  on  each  subcycle.  It  is  therefore  possible  to  create  a 
mixed  cell  after  one  or  more  subcycles  of  INFACE  have  been 
completed.  In  that  case  the  flux  terms  (SAMMP,  SAMPY)  for 
that. cell  have  not  been  accumulated  for  the  completed  sub-  . 
cycles  and  need  to  be  computed  before  adding  in  the  flux 
terms  for  the  current  and  subsequent  subcycles.  Therefore, 
when  INFACE  has  completed  at  least  one  subcycle,  NEWMIX  calls 
FLUX  after  setting  up  the  storage  for  the  new  mixed  cell. 
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4.7  CREATING  PURE  CELLS 


4.7.1  Defining  Material  that  Remains  in  the  Cell 

After  the  fractional  areas  associated  with  each 
material  package  have  been  computed,  INFACE  searches  for 
cells  that  were  interface  cells  on  the  previous  cycle  but 
no  longer  are  cut  by  an  interface.  Such  a  cell  has  become 
pure  and  its  flag,  MFLAG,  is  made  negative  to  signify  that 
fact  until  the  transport  has  been  completed  in  PH2.  To 
determine  what  material  package  a  new  pure  cell  K  belongs 
to,  INFACE  first  looks  for  a  neighbor  that  is  pure  and 
assumes  cell  K  will  belong  to  the  same  package  as  one  of 
its  pure  neighbors,  e.g.,  in  the  following  diagram, 


the  cell  below  is  pure  (MFLAG(KB)  =  2) ,  therefore  cell  K  has 
become  a  pure  cell  be  longing  to  package  2.  In  case  all  four 
neighbors  are  mixed,  the  cells  KB  and  KL  are  examined  as 
illustrated  in  the  following  example. 
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Void 
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KB 
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X 


Material  Package  1 


X 


\ 


Material  Package  2 


Since  cell  K  is  not  crossed  by  an  interface  it  follows  that 
if  the  fractional  area  of  the  top  of  cell  KB  for  package  2 
is  the  total  cell  boundary  area  then  cell  K  must  belong  to 
package  2.  The  cell  on  the  left  of  cell  K  is  used  analogously 
if  cell  K  is  in  the  bottom  row  of  the  grid. 

4*7*2  Evacuating  Materials 

To  signify  that  a  material • interface  has  passed  out 
of  cell  K,  the  density  of  that  material  is  set  to  zero.  For 
example,  consider  the  figure  above  assuming  there  are  two 
material  packages  in  the  grid;  the  densities  of  cell  K  are 
redefined  until  the  end  of  the  cycle  as  follows: 

M  =  MFLAG(K)  -  100 

RHO (1 ,M) '=  0. 

RHO (2 ,M)  f  0. 
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These  zero  densities  are  later  used  as  a  signal  to  adjust  the 
fluxes  of  material  1  so  that  the  material  will  be  exactly  eva¬ 
cuated  from  cell  K.  This  adjustment  of  the  fluxes  usually  occurs 
twice  in  INFACE  by  calling  subroutine  ADJFLX,  once  during  the 
first  subcycle  (for  cells  that  become  pure  during  the  first  sub¬ 
cycle)  and  once  after  all  subcycles  are  completed  (for  cells 
that  become  pure  after  the  first  subcycle) .  For  the  first 
case  the  flux  terms  for  mixed  cells  (SAMMP,  SAMPY)  are  saved 
from  the  previous  cycles  to  indicate  the  direction  of  the 
flow  and  the  direction  of  the  evaluation.  In  the  following 
example,  t  represents  time,  and  ML  and  MB  are  locations 
associated  with  the  mixed  cells  on  the  left  of  and  below  cell 
K,  respectively: 

t 

2 


'A 

X 

Ax 

_  X 

Material 

Package 

N 

KL 

s 

X 

X  V 

X 

K 

.  X 

X  \ 

X  ' 
s 

\  X 

* 

Void 

"  KB 

X 

Material 

Package 

S. 

t  flux  terms 
1 

SAMMP (1 ,M)  >  0. 
SAMPY (1,M)  >  0. 
SAMMP (1, ML)  =  0. 
SAMPY (1, MB)  =  0. 


At  t  ,  the  t  flux  terms  indicate  that  any  mass  of 

2  i 

package  1  still  in  cell  K  should  be  evaluated  out  the  right 
and  top  boundaries,  not  the  left  or  bottom.  These  evacua¬ 
tion  procedures  are  used  for  a  cell  that  has  lost  one  inter¬ 
face  but  is  still  cut  by  another  interface  (so  is  still 
mixed)  as  well  as  for  cells  that  have  become  pure. 


56 


V.  THE  PRESSURE  AND  DENSITY  ITERATION  FOR  MIXED. CELLS 
5.1  GENERAL  METHOD 

For  ordinary  one -material  cells,  it  is  possible  to’ 
compute  the  pressure  directly  from  the  material  density  and 
specific  internal  energy.  For  cells  containing  more  than 
one  material,  it  'is  necessary  to  make  additional  assumptions 
in  order  to  determine  cell  pressure.  In  the  present  method, 
the  cell  pressure  for  mixed  cells  is  determined  by  an  intera- 
tion  procedure.  Specifically,  the  densities  of  the  various 
materials  within  the  cell  are  varied,  subject  to  the  restric¬ 
tion  that  the  cell  be  exactly  filled  by  the  masses  therein, 
until  the  individual  pressures  converge  to  a  common  value, 
taken  to  be  the  cell  pressure.  This  process  gives,  in  addi¬ 
tion  to  the  cell  pressure,  the  densities  of  the  individual  • 
materials  for  subsequent  use  in  the  Phase  1  energy  partition 
and  in  the  Phase  2  transport  calculation. 

The  procedure  for  determining  the  pressure  and  the 
densities  in  a  mixed. cell  is  actually  divided  into  two  parts, 
a  pre-iteration  calculation  to  assure  that  the  cell  is  filled 
exactly  and  an  iteration  calculation  to  converge  on  the  de¬ 
sired  P's  and  p  Ts .  - 
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5.2  PRE-ITERATION  CALCULATION 


In  general,  the  masses  and  densities  which  exist  at 
the  beginning  of  the  calculation  do  not  exactly  fill  the 
cell  since  cell  masses  were  changed  in  the  Phase  2  calcula¬ 
tion  of  the  previous  cycle.  The  pre - iteration  calculation 
fills  the  cell  exactly  by  taking  account  of  the  constant- 
energy  compressibilities  of  each  of  the  various  materials. 
Stepwise,  this  is  done  by: 


1. 


Calling  EQST  to  compute  and  C| 

3P^  /  3  p^  from  the  old  values 

o f  and  (assumed  constant  through¬ 

out  the  calculation)  for  each  material  in  the 
cell. 


2..  Normalizing  (see  proof  below)  the . dens i ties 
to  fill  the  cell  exactly 

A V i  =  -AP/h? 


where 


AP 


VOL 


_  (VOL-VCELL) 


i 
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5.3  ITERATION  CALCULATION 


•  Given  that  p.  =  f ■  (v.  ,  E.)  and  C?  -  (sP- /9p-  )F . 
(calculated  by  EQST)  for  materials  i  =  1 ,  .  . .  M,  it  is 
desired  to  find  cell  pressure,  P,  by  varying  the  (=  1/p^ 

until  the  P ^ s  are  equal  within  some  specified  accuracy, 
subject  to  the  -restriction  that  the  cell  remains  exactly 
filled  and  that  the  are  held  constant. 

The  equation  of  volume  conservation  is 

51  nu  AVt  =  0 


and  the  equation  that  causes  material  i  to  undergo  a 
change  in  specific  volume  V^,  such  that  its  pressure  change 
from  its  current  value  P^  =  f^V^,  )  to  a  value  P^+^ , 

common  to  all  materials  in  the  cell  at  the  end  of  the 
iteration  cycle,  is  . 


=  (P«+1 


where 


These  equations,  for  i  =  1,  2,  M  can  be  regarded  as 

M+l  equations  for  M+l  unknown  quantities,  AV^  and  P^+^. 
Other  quantities  in  the  equations  are  either  known  constants 
(m-)  or  are  updated  each  cycle  of  the  iteration 
|p?  and  (3V./3P^)^  I  and  are  taken  to  be  constants  while 
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The  solutions 


solving  these  equations  for  AV^  and  PN+1. 
are  1 


pN+l 


Ewi 


and 


W.  = 

i 


Vhi 


That  these  equations  are  solutions  to  the  given  equations 
can  be  verified  by  direct  substitution.- 


The  iteration  would  be  exact  (single  cycle  convergence) 
if  the  input  coefficients  ( 3 P ^ / 3 V ^ ^  were  constants,  since 
the  solution  does  not  involve  additional  approximations . 
However,  since  this  is  not  the  case,. the  coefficients  and 
the  pressures  are  necessarily  recomputed  each  step  of  the 
iteration  by  calling  subroutine  EQST  with  the  V-  generated 
during  the  previous  step. 


Stepwise,  the  iteraction  calculation  is  as  follows: 

1.  Call  EQST  to  find  Pi  =  fi(vi,  E±) 
and  (3V./3P.)e_ 
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2.  Compute  new  P^+^ 

3.  Compute  new  AV^  from  above  formulas . 

4.  Update  by  adding  the  above  AV^  to 

the  old  specific  volume. 

5.  Return  to  (1)  using  the  new  values  for  . 

The  iteration  is  complete  when  the  P^’s  all  equal 
P  to  some  preset  accuracy.  Usually,  convergence  is  obtained 
in  a  maximum  of  two  or  three  steps  if  the  convergence  cri¬ 
terion  |P^  -  P|  <  10  3  P  is  used. 
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VI.  THE  PRESSURE  BOUNDARY  CONDITION 


The  BRLSC  code  has  the  option  of  applying  a  time  and 
position  dependent  pressure  boundary  condition. 

6.1  GENERAL  METHOD 

When  a  pressure  boundary  condition  is  to  be  applied, 
subroutine  APLYPR  is  called  at  the  beginning  of  CDT.  This 
subroutine  identifies  the  section  of  the  moving  surface 
along  which  the  pressure  force  is  to  be  applied  and  determines 
which  cells  are  cut  by  that  portion  of  surface.  Those  cells 
are  flagged  and  the  coordinates  of  a  point  in  the  cell  which 
lies  on  the  surface  are  determined.  These  coordinates  ,  along 
with  the  time,  are  used  as  input  to  subroutine  PRESBC*  which 
determines  the  pressure,  density  and  internal  energy  of  the 
driving  material  at  that 'position,  and  time.  The  pressure, 
density  and  internal  energy  are  then  returned  to  APLYPR 
which  iterates  on  the  density  of  the  material  already. in  the 
cell  until  the  pressure  of  that  material  is  equal  to  the  pres¬ 
sure  which  was  returned  from  PRESBC.  From  the  mass  and  new 
density  of  the  material  in  the  cell,  the  volume  which  the 
material  occupies  is  determined  and  the  volume  of  the  void 
remaining  in  the  cell  computed.- .  The  volume  of  the  void 
is  then  filled  with  mass  with  the  density  and  internal 
energy  which  was  returned  from  PRESBC;  The  total  cell  mass 
and  internal  energy  and  the  total  theoretical  energy  in  the 
problem  are  updated  to  reflect  this  extra  mass  and  energy 
and  the  cell  pressure  is  set  equal  to  the  desired  driving 
pressure.  When  all  cells  cut  by  the  pressure- loaded  section 
of  the  free  surface  have  been  flagged,  the  pressures  set  and 
the  masses  and  energies  updated,  APLYPR  returns  to  CDT. 

The  user  must  supply  subroutine  PRESBC  if  a  pressure 

boundary  condition  is  to  be  employed. 
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Any  cells  which  have  had  their  pressures  set  in  APLYPR  are 
ignored  during  the  pressure  iteration  section  of  CDT  but 
are  included  when  determining  the  time  step. 

After  the  pressures  have  been  defined  along  the 
desired  section  of  the  material  boundary,  they  are  used  to 
accelerate  the  material.  This  is  done  in  PHI  by  calling  sub¬ 
routine  FAKMAT  for  all  cells  which  have  been  flagged  by 
APLYPR.  FAKMAT  looks  at  the  four  cells  which  surround  the 
cell  being  processed.  If  a  neighboring  cell  is  flagged  as 
being  a  void  cell,  the  interface  between  that  cell  and  the 
cell  (K)  being  processed  by  PHI  is  considered  to  be  a  trans- 
mittive  boundary.  The  interface  pressure  (which  PHI  had  set 
equal  to  zero),  is  reset  equal  to  the  pressure  of  cell  K  and 
the  work  done  at  the  interface  is  added  to  the  total 
theoretical  energy  of  the  problem.  PHI  then  utilizes  the 
new  values  for  interface  pressures  to  update  the  cell’s 
velocity  and  internal  energy.  The  energy  is  then  partitioned 
between  the  original  material  in  the  cell  and  the  artifi¬ 
cially  added  driving  material  which  ivas  added  by  APLYPR. 

The  last  step  in  applying  the  pressure  boundary 
condition  is  to  call  RMVMAT  at  the  end  of  PHI.  This  sub¬ 
routine  removes  the  mass  and  energy  of  the  material  which 
was  added  by  APLYPR,  redefines  the  total  cell  mass  and 
internal  energy,  and  updates  the  total  theoretical  energy 
to  reflect  the  reduction  in  total  energy.  The  flags  that 
were  set  by  APLYPR  are  unset  and  RMVMAT  returns  control 
to  PHI.  At  this  point,  the  grid  again  contains  only  one' 
material  but  its  total  energy  and  momentum  reflect  the 
work  and  impulse  delivered  it  by  the  pressure  boundary  con¬ 
ditions  . 

Subroutine  RMVMAT  is  also  called  from  one  other  place 
in  the  code.  Since  restart  tape  dumps  arc  made  in  EDIT,  and 
EDIT  is  called  after  CDT  and  before  Pill,  the  cell  masses  and 
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internal  energies  written  on  the  tape  include  the  artifically 
added  mass  and  energy  added  by  APLYPR.  It  is  therefore 
necessary  to  call  RMVWAT  after  INPUT  and  before  CDT  any  time 
a  problem  is  being  restarted  in  order  to  restore  the  grid.  to. 
the  pre-CDT  condition  it  had  before  the  tape' dump  was  made. 

The  following  section  shows  step  by  step  how  the 
pressure  boundary  condition  is  applied. 


STEP-BY-STEP  APPLICATION  OF  THE  PRESSURE  BOUNDARY 
CONDITION 


Free  Surface 
Interface 


6.2.1  Subroutine  APLYPR 

Subroutine  APLYPR,  called  at  the  beginning  of  CDT, 
follows  the  interface  tracers  which  define  the  portion  of 
the  material  boundary  along  which  the  pressure  force  is  to 
be  -applied  and  determines  which  cells  it  intercepts  (in  the 
above  sketch,  the  interface  enters-cell  K  at  Y^) 

and  leaves  at  (X  .  ,  Y  .  )  .  Then,  for  each  cell  K: 


(1)  Compute 


(3)  Iterate  on  p^  (the  density  of  liner  material 
in  cell  K)  until 

P£  =  EP  =  PX 

where  E£.  =  specific  internal  energy  of  the  liner 
material  and  P£  is  the  pressure. 

(4)  Use  the  updated  value  for  to  compute  the 

volume  of  X  which  will  be  filled  with  explosive 

Volx  =  VCELL  -  m£/Pjl 

(5)  Fill  the  cell  with  mass  (with  values  defined  by 
PRESBC)  and  update  the  total  theoretical  energy 
to  reflect  this  added  mass  and  energy 

Am  =  p^  •  Vol^ 

ETH  =  ETH  +  AM ^  |  (U(K)2  +  V(K)2) 

AIX(K)  =  (AMX(K)  •  AIX(K)  +  Am  •  EX)/(AMX(K)  +  Am) 
AMX(K)  =  AMX(K)  +  Am 
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(6)  Define  pressure  in  cell 

P(K)  =  Px 

6.2.2  Pressure  Calculations  in  CDT 

When  control  returns  from  APLYPR  to  CDT,  the  pressures 
in  all  cells  along  which  the  pressure  force  is  to  be  applied 
have  been  defined.  Also,  the  masses ,  and  energies  in  these 
cells  have  been  updated  to  simulate  the  existence  of  the  ex¬ 
plosive  material  in  the  cells.  CDT  skips  these  cells  when 
computing  cell  pressures  since  their  pressures  have  already 
been  computed. 

6.2.3  Special ' Treatment  for  Pressure  Loaded  Cells  in  PHI 

(1)  Special  consideration  must  be  given  to  any  cell 
interface  which  exists  between  a  cell  which  has  had  its  pres¬ 
sure  defined  by  APLYPR  and  a  void  cell.  It  must  be  remem¬ 
bered  that  if  the  driving  material  were  actually  in  the 
problem  instead  of  being  defined  by  a  function  (PRESBC) , 
cells  bordering  the  interface  cells  that  have  had  their 
pressures  set  would  not  be  void  but  would  contain  the 
driving  material.  It  then  becomes  necessary  to  treat  these 
special  cell  interfaces  as  being  interfaces  between  two  cells 
containing  material  instead  of  as  an  interface,  between  a 
cell  containing  material  and  a  void  cell.  This  is  done  by 
subroutine  FA KM AT  which  identifies  these  special  cell  inter¬ 
faces  and  sets  their  values  of  velocities  and  pressures  so 
that  these  cells  can  be  treated  as  if  they  were  trans- 
mittivc  boundaries.  Referring  to  the  previous  sketch,  con- 
side' r  the  interface  between  cell  K  and  void  cell  KR. 

During  the  normal  execution  of  Phase  1,  the  following  inter¬ 
face  values  would  be  computed  as  (see  Section  2.2,2): 
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and  energies  in  cell  K  are 
updated  in  PHI,  subroutine  FAXMAT  would  be  called.  FAKMAT 
would  redefine  the  interface  pressure: 

Pr  =  P(K) 

and  update  the  total  theoretical  energy  in  the  problem  by 

ETH  =  ETH  +  Ar  •  Ur  *  Pr  •  At  . 
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Then,  before  the  velocities 


The  interface  value  of  pressure  between  cell  K  and  void 
cell  KA  would  also  be  modified  in  a  similar  manner. 

(2)  Once  the  total  change  in  internal  energy  for 
the  entire  cell  has  been  computed  (AAIX)  it  is  necessary 
to  increment  the  internal  energy  in  the  liner  (ASIE  ),  This 
is  done  by  applying  the  following  equation 


ASIE^  = 


AMX(K)  *  AAIX 


Mi  (AMX(K)  -  M£) 


(piC£  PxYxp(X) 


where  and  yx  were  determined  in  PRESBC  and  is 

the  sound  speed  in  the  liner. 


(3)  After  the  applied  pressure  force  has  been  used 
to  update  the  velocities  and  energies  in  tiie  interface  cells, 
the  mass  and  associated  energy  which  was  added  by  APLYPR 
must  be  removed.  This  is  done  by  calling  .RMVMAT  at  the  end 
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of  Pill.  RMVMAT  redefines  the  total  theoretical  energy, 

ETH  =  ETil  +  E£  -  AMX(k)  •  AIX(k)  +  |[m4  -  AMX(k)J  [  U(K)  :+V(K) 

and  redefines  the  total  cell  mass  and  energy 

AMX(K)  -  M 
AIX(k)  -  E£ 

At  the  end  of  PHI,  the  momentum  and  energy  of  the 
liner  material  has  been  updated  to  reflect  the  impulse  and 
work  delivered  to  it  by  the  applied  pressure  boundary  con¬ 
dition. 


70 


VII'.  SAMPLE  CALCULATIONS 


In  this  section  of  the  report,  selected  results  from 
three  calculations  are  given.  The  three  calculations  were 
chosen  to  demonstrate  the  three  modes  in  which  problems  can 
be  run  with  BRLSC  (see-  Section  8. 3. 1.1).  For  a  more  de¬ 
tailed  description  of  several  problems  run  by  BRLSC,  see 
Ref.  [7]. 

7.1  THE  BRL  PRECISION  SHAPED  CHARGE 

* 

This  problem  was  run  in  mode  1  using  the  Shaped 
Charge  reported  in  Ref.  [4] .  The  explosive  was  octol  and 
the  liner  was  copper  with  a  thickness  of  .081  inches  and 
a  half  angle  of  21  degrees.  The  liner  was  set  up  in  an 
Eulerian  grid  with  square  cells  (DX  =  DY  =  .145  cm). 

Strength  effects  were  not  included  in  the  calculation 
since  earlier  calculations  had  shown  that  including  strength 
in  shaped  charge  calculations  has  a  negligible  effect  on 
the  resulting  jet. 

Figures  4  and  5  show  the  conf igurations  at  various 
times.  For  further  details  about  this  problem,  see  Ref.  [7] 

7 . 2  AN  ARTIFICIALLY  DRIVEN  SHAPED  CHARGE  LINER 

*  * 

This  problem  was  run  to  demonstrate  the  results 
obtained  by  running  a  problem  in  mode  2  (using  a  time  and 
position  dependent  pressure  force  along  the  outer  surface 
of  the  liner).  The  liner,  a  shortened  version  of  the  105  mm 
shaped  charge  liner^,  was  copper,  .106  inches  thick,  and 
had  a  half  angle  of  21  degrees.  The  grid  contained  square 

X 

This  problem  had  5600  cells  in  the  Eulerian  grid  and  took 
102  minutes  to  run  471  cycles  (55.0  usee)  on  a  Univac  1108 

■k  ± 

This  problem  had  2100  cells  in  the  Eulerian  grid  and  took 
10.5  minutes  to  run  196  cycles  (10.0  usee)  on  a  CDC  6600. 
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Fig.  4--BRL  precision  shaped  charge. 
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cells  (DX  =  DY  =  .12  cm)  and  the  effects  of  strength  were 
included. 

Figure  6  sho\\?s  the  liner  configuration  at  10  ysec. 

The  shaded  regions  shows  cells  which  have  not  failed.  For 
additional  details,  see  Ref.  [7]. 

7.3  THE  COLLAPSE  OF  A  FREE  FLYING  LINER 

This  problem  was  run  during  the  early  stages  of 
developing  BRLSC  and  is  included  as  an  example  of  mode  3. 

The  liner  was  set  up  with  an  initial,  velocity  of  3.3  *  10s 
cm/sec  directed  inward  in  the  direction  of  the  liner.  The 
liner  was  copper,  had  a  thickness  of  0.27  cm,  and  had  a  half 
angle  of  35.75  degrees.  The  grid  was  made  up  of  square  cells 
with  (DX  =  DY  =  .05  cm).-  The  effects  of  strength  were  not 
included. 

Figure  7  shows  the  configurations  at  various  times. 

For  additional  details,  see  Ref.  [6]. 
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t  =  0 


t  =  0,6  usee 


Fig.  7- -Example  of  a  Mode  3  calculation 
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VIII.  USERS  GUIDE  TO  THE  BRLSC  PROGRAM 


This  section  of  the  report  constitutes  a  users  manual. 

It  contains  the  information  necessary  to  generate  and  run  shaped 
charge  problems  with  the  BRLSC  code  and  is  divided  into  five 
parts.  The  first  part  gives  a  short  description  of  some  of 
the  major  problems  encountered  while  writing  BRLSC  and  what 
was  done  to  overcome  or  minimize  these  problems.  The  second 
part  contains  a  brief  description  of  every  subroutine  in  the 
BRLSC  code.  The  third  part  gives  some  of  the  code's  restric¬ 
tions,  the  options  available,  and  explains  how  to  set  up  the 
input  deck  for  generating  and  running  shaped  charge  problems. 

In  Section  8.4,  the  procedure  for  restarting  problems  is  given 
and  Section  8.5  is  a  dictionary  of  the  more  important  variables 
used  in  BRLSC. 

8.1  INTRODUCTORY  REMARKS 

BRLSC  is  basically  the  HELP  ^  code  which  has  been 
modified  to  allow  it  to  handle  shaped  charge  problems.  While 
writing  BRLSC,  it  became  apparent  that  there  were  several  main 
problems  to  overcome.  A  brief  discussion  of  the  problems  and 
the  modifications  made  to  correct  them  follows. 

8.1.1  Problem  Generation 

It  was  necessary  to  rewrite  subroutine  SETUP  to  allow 
shaped  charge  problems  to  be  easily  generated,  giving  the  user 
a  wide  choice  of  possible  materials,  configurations,  zoning, 
etc. 
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8.1.2  Thin  Liners  and  Free  Surface  Cells 


Since  shaped  charge  problems  involve  thin  liners  moving 
diagonally  through  the  computational  grid,  the  majority  of  the 
cells  which  contain  liner  material  are  mixed  cells,  either  con¬ 
taining  liner  material  and  high  explosive  or  liner  material 
and  void.  Two  main  problems  became  apparent  when  a  thin  liner 
was  moved  through  the  grid.  First,  free  surface  or  underdense 
cells  tended  to  be  over-accelerated  in  PHI.  This  problem  has 
been  solved  by  rewriting  the  PHI  difference  equations  so  that 
interface  pressures  are  computed  as  an  inverse  density  weighted 
average  of  the  cell  pressures.  The  second  problem  involved  the 
realistic  definition  of  the  material  density  in  free  surface 
cells  used  in  computing  mass  fluxes.  This  problem  was  overcome 
by  computing  volume  fluxes  as  well  as  mass  fluxes  and  updating 
the  density  of  the  material  in  free  surface  cells  at  the  end 
of  PH2.  Also,  the  material  densities  are  updated  in  PHI  by 
computing  the  change  in  volume  due  to  the  surrounding  velocity 
field. 

8.1.3  Computational  Stability 

In  the  initial  solutions  of  shaped  charge  problems,  it 
was  observed  that  instabilities  occurred  in  regions  of  the 
grid  with  low  velocities.  These  instabilities  often  lead  to 
unsatisfactory  results.  The  stability  of  a  calculation  depends 
in  part,  on  the  effective  viscosity  which  is  the  result  of 
converting  kinetic  energy  to  internal  energy  during  mass  trans¬ 
port  (PH2) .  If  velocities  are  small,  there  is  little  mass  flux 
between  cells  and,  as  a  result,  little  effective  viscosity. 
Thes.e  instabilities  were  considerably  reduced  by  computing  a 
pseudo-viscous  pressure  term  (i.e.,  artificial  viscosity) 
which  is  included  in  the  interface  pressures  computed  in  PHI. 
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8.1.4  Negative  Internal  Energies 


Another  major  problem  which  arose  during  the  early  stages 
of  the  development  of  BRLSC  involved  negative  internal  energy. 
Negative  internal  energies  were  generated  in  one  of  three  ways. 
First,  when  cells  were  over- accelerated  in  PHI,  the  increase 
in  kinetic  energy  in  the  underdense  "kicked''  cell  was  offset 
by  a  decrease  in  internal  energy  in  the  cell(s)  which  did  the 
"kicking".  Secondly,  if  the  mass  fluxes  out  of  a  cell  were 
greater  than  the  mass  in  the  cell  (regardless  of  the  mass 
fluxes  into  the  cell),  negative  internal  energies  could  result. 
Finally  the  internal  energy  of  one  material  in  a  mixed  cell 
could  become  negative  when  its  momentum  and  energy  were  updated, 
even  if  the  total  cell  internal  energy. was  positive.  This  can 
occur  in  PHI  or  PH2 .  The  first  of  these  problems  t\ras  minimized 
by  the  inverse  density  weighting  scheme  for  computing  inter¬ 
face  pressures  in  PHI.  The  second  was  cured  by  placing  limi¬ 
tations  on  the  mass-  fluxes  out  of  cells.  The  third  was  solved 
by  modifying  the  manner  in  which  the  internal  energies  of 
materials  in  mixed  cells  are  incremented  in  PHI  and  by  insuring 
that  all  internal  energies  of  materials  in  mixed  cells  have 
the  same  sign  as  the  total  cell  energy  at  the  end  of  PHI  and 
PH2 . 

8.1.5  Overheating  of  Cells 

Certain  cells  in  regions  of  very  high  ve loci ty  gradients 
tend  to  become  overheated,  i.e.,  their  calculated  internal  energy 
is  unrealistically  high.  This  arises  when  the  mass  entering  a 
cell  is  moving  at  a  velocity  extremely  different  from  the 
mass  already  in  the  cell,  so  that  too  much  kinetic  energy  is 
thermalized.  To  prevent  this  overheating  from  unrealistically 
vaporizing  material,  only  the  condensed  form  of  the  Tillotson 
equation  of  state  is  used. 
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8.1.6  Explosives 


In  order  to  run  shaped  charge  problems,  it  was  necessary 
to  modify  the  code  to  handle  explosives.  This  was  done  by 
adding  subroutine  EOSXPL  which  is  called  by  EQST  when  computing 
pressures  in  cells  containing  the  explosive.  EOSXPL  determines 
the  location  of  the  detonation  front  and  whether  the  cell  in 
question  is  behind,  intercepted  by,  or  ahead  of  the  detonation 
front,  modifies  the  equation  of  state  for  the  explosive 
accordingly,  and  computes  the  pressure. 

8.1.7  The  Pressure  Boundary  Condition 

Shaped  charge  problems  can  be  run  with  BRLSC  by  applying 
time  and  position  dependent  pressure  force  along  the  outer 
surface  of  the  liner.  Subroutine  PRESBC  must  be  supplied  the 
necessary  information  to  compute  the  explosive's  pressure, 
density,  and  internal  energy  for  any  point  along  the  liner- 
explosive-  interface  for  any  time.  These  values  can  be  deter¬ 
mined  by  a  crudely  resolved  BRLSC  calculation,  a  Lagrangian 
calculation,  or  from  analytical  approximations.  BRLSC  then 
computes  the  resulting  liner  collapse  and  jet  formation.  Since 
there  is  no  actual  explosive  material  in  the  grid,  the  cost  of 
running  a  problem  involving  a  well  resolved  liner  is  greatly 
reduced.  Also,  by  omitting  explosive  material  from  the  problem, 
an  effective  slip  surface  exists  between  the  would-be  explosive- 
liner  interface. 
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8.2  SUBROUTINE  DESCRIPTIONS 


The  following  section  contains  a  brief  description 
of  every  subroutine  used  by  the  BRLSC  code.  Also  included 
is  a  general  flow  diag-ram  (Fig.  9)  showing  the  inter-relation¬ 
ship  between  subroutines,  and  Table  3  lists  all  of  the  sub¬ 
routines,  which  subroutines  they  call,  and  which  subroutines 
they  are  called  from. 

ADJFLX:  Subroutine  ADJFLX  calculates  the  mass  and 
volume  fluxes  out  of  mixed  cells  which  have  lost-  a  material 
interface  and  should  be  evaluacted  of  that  material.  It  also 
initializes  the  mixed  cell  fluxes  when  called  during  the  first 
pass  through  INFACE. 

ADJSIE :  Subroutine  ADJSIE  makes  sure  that  all  the 
specific  internal  energies  assigned  to  materials  in  a  mixed 
cell  have  the  same  sign  as  the  total  specific  internal  energy 
for  the  cell.  It  also  insures  that  the  sum  of  the  cellTs 
material  masses  and  energies  equal  the  total  cell  mass  and 
energy. 

APLYPR :  Subroutine  APLYPR  can  be  used  to  apply  a  time 

and  position  dependent  pressure  boundary  condition  to  a  section 
of  a  material  boundary.  It. identifies  the  cells  to  which  the 
pressure  is  to  be  applied,  calls  PRESBC  to  determine  the 
pressure,  density  and  specific  internal  energy  that  the  arti¬ 
ficially  applied-  material  has  at  that  point,  adjusts  the 
density  of  the  material  already  in  the  cell  to  give  the 
desired  pressure,  and  fills  the  rest  of  the  cell  with  the 
driving  material. 
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Subroutines  which  have 
calls  to  LRliott 


Fig,  9- -General  flow  diagram  of  BRLSC  program. 
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Subroutine 


EDITOR 


EOSXPL 


ERROR 


FAKMAT 


:  TABLE  3 

SUBROUTINE  REFERENCES 


EQST,  NEiVMIX,  PRESBC 

CDT,  EDIT,  EDITOR, 
INFACE,  INPUT,  PHI,  PH2 
RMVMAT 

ERROR,  NEIVMIX,  NEWRHO 


APLYPR,  EQST,  ERROR 
ERROR,  MAP 


EOSXPL 


EDIT 


Is  Called  From 


INFACE 
PHI,  PH2 


INFACE 

INPUT 


BRLSC 


BRLSC,  ERROR 

BRLSC,  SETUP,  SETUPS 

EQST,  SETUP 

APLYPR,  CDT 

CALFRC ,  CDT,  EDIT,  • 
INFACE,  INPUT,  NEIVMIX, 
NEWRHO,  PH2 ,  SETUP 


INFACE,  NEiVMIX 


INFACE 


LOCI  J 


ADJFLX,  CALFRC,  ERROR, 
FLUX,  MOVTCR 

CARDS,  ERROR,  SETUP 


BRLSC 


BRLSC 


SETUP 2,  SETUP 4 ,  SETUP5 


EDIT 


MOVTCR 


INFACE 


NEiVMIX 


ERROR,  FLUX 


APLYPR,  CALFRC,  SETUP2 


Table  3 

Subroutine  References  (continued) 


Subroutine 

Calls 

Is  Called  From 

NEWRHO 

ERROR 

CALFRC 

PHI 

ADJSIE,  FAKMAT,  PH3 , 
PRSVIS,  RMVMAT 

BRLSC 

PH  2 

ADJSIE,  ERROR,  ZROMAS 

BRLSC 

PH3 

STRNG 

PHI 

PRESBC 

APLYPR 

PRSVIS 

PHI 

RMVMAT 

BRLSC,  PHI 

SETUP 

EDITOR,  ERROR,  SETUP1 
SETUP4 ,  SETUP5 ,  SETXPL* 

I  INPUT 

■ 

SETUP  1 

SETUP2 

SETUP 

SETUP2 

LOCI J ,  NEWMIX,  SETUP3 

SETUP1 

SETUP3 

SETUP2 

SETUP4 

LOCI  J 

SETUP 

SETUP 5 

EDITOR,  LOCIJ 

SETUP 

STRNG 

PH3 

Z ROMAS 

PH2 

SETXPL  is  an  entry  point  in  EOSXPL 
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BLOCK :  A  block  data  element  that  defines  normal 

density  and  the  bulk  sound  speed  of  nineteen  materials. 

BRLSC :  The  overall  cycling  of  the  calculation  is 
controlled  by  BRLSC,  as  shown  in  Fig.  9.  It  is  the  main 
routine  of  the  code  and  calls  INPUT,  RMVMAT,  CDT,  EDIT,  PHI, 
INFACC  and  PH2 ,  in  that  order.  If  it  is  necessary  to  see 
the  effects  of  each  phase  of  the  calculational  cycle,  BRLSC. 
will  respond  to  the  input  parameter  INTER  and  call  EDIT  on 
print  cycles  after  PHI  as  well  as  after  CDT.  BRLSC  also 
calls  EDITCR  and  EXIT  on  the  normal  cessation  of  the 
calculation . 

CALFRC :  In  subroutine  CALFRC,  the  tracer  particles 

that  circumscribe  each  material  package  are  followed  and 
the  intercepts  with  cell  boundaries  are  determined.  The 
location  of  each  intercept  is  then  used  to  determine  the 
fractional  cell  areas  for  each  material  in  each  mixed  cell. 
If  a  previously  pure  cell  is  intercepted,  it  is  flagged 
mixed. 


CARDS :  Most  of  the  input  parameters  stored  in  blank 
common  are  read  by  subroutine  CARDS.  The  format  of  these 
input  cards  is  described  in  Section  8.3.2. 

CDT :  A  principal  function  of  this  routine  is  to  compute 
a  time  step  which  ensures  stability  of  the  finite  difference 
equations.  This  is  done  by  finding  the  cell  with  the  minimum 
of 

(a)  AXAY/ (AX  1 V |  +  AY | U | ) 

(b)  AX/C 

(c)  AY/C 

where  AX  is  the  radial  dimension  of  the  cell  (cm) ,  AY  is 
the  axial  dimension  of  the  cell  (cm),  |U|  is  the  absolute 
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value  of  the  radial  velocity  of  the  cell  (cm/scc) ,  V  is  the 
absolute  value  of  the  axial  velocity  of  the  cell  (cm/sec)  and 
C  is  the  local  sound  speed  (cm/sec) .  For  a  polytropic  gas 
the  sound  speed  is  computed  as  /yp/p  and  for  other  materials 
by  the  appropriate  relation  C  =  +  B  /P  where  P  is  the 

pressure  in  the  cell  and  the  coefficient,  IT,  is  an  input  para¬ 
meter.  C  is  defined  in  BLOCK  for  nineteen  materials  as 

/A/p  ,  where  A  is  a  coefficient  in  the  solid  equation  of 
o 

state.  In  a  mixed  cell  the  sound  speed  is  given  by  a  mass 
weighted  average  of  the  sound  speeds  of  the  materials  in  the 
cell.  Each  cycle  ' CDT  prints  the  column  and  row  (I,J)  of 
the  cell  controlling  the  time  step  as  well  as  the  maximum  sound 
speed  in  the  grid  (MAXC) ,  the  maximum  velocity  in  the  grid 
(jMAXUV)  and  the  velocity  and  pressure  cutoff  values. 

Another  function  of  CDT  is  to  equilibrate  the  pressures 
of  materials  in  mixed  cells  using  an  iteration  scheme  that 
adjusts  the  material  densities.  A  detailed  discussion  of  this 
iteration  method  is  found  in  Section  5.  The  pressures  for 
pure  cells  are  also  updated  in  CDT  by  a  call  to  the  equation 
of  state  subroutine,  EQST,  with  the  density,  energy  and 
material  code  number  of  the  cell  (RHOW,  ENERGY,  N)  passed 
through  blank  common. 

When  a  pressure  boundary  condition  is  being  applied  to 
a  material  boundary,  CDT  calls  subroutine  APLYPR  which  com¬ 
putes  the  necessary  densities  and  pressures  in  the  cells  that 
the  surface  intercepts.  These  cells  are  bypassed  in  the 
pressure  iteration  section  of  CDT. 

EDIT :  The  periodic  printing  and  writing  on  the  re¬ 

start  tape  are  executed  by  subroutine  EDIT.  The  frequency 
of  printing  and  tape  dumps  is  controlled  by  input  parameters 
described  in  Section  8.3.3. 
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On  every  print  cycle  EDIT  prints  the  mass,  total  energy, 
internal  energy,  kinetic  energy,  axial  and  radial  momentum, 
and  plastic  work  for  each  material  package  as  well  as  for  the 
entire  grid.  The  changes  in  energy  due  to. evaporation  and 
losses  out  boundaries  are  also  accounted  for.  The  coordinates 
of  the  material  tracers  circumscribing  each  package  are  printed 
in  cell  units.  Summary  maps  of  the  compression,  pressure, 
velocities  and  internal  energies  are  printed  if  the  input 
parameter,  MAPS,  is  non-zero.  And  finally,  the  pressure, 
velocities,  internal  energy,  compression  and  stress  deviators 
for  all  cells  in  the  active  grid  are  displayed  on  "long"  EDIT 
prints  and  for  each  cell  in  the  axis  column  on  "short”  EDIT 
print. 

Another  important  function  of  EDIT  is  to  compute  the 
relative  error  in  the  total  energy  of  the  grid  when  summed 
over  all  the  cells,  and  to  check  that  this  error  does  not 
exceed  a  limit  specified  by  the  input  variable  DMIN. 

EDIT  also  senses  when  execution  should  stop  and  sets 
the  exit  flag,  WFLAGL. 

EDITCR :  The  coordinates  of  the  interior  tracers  (XP , 

YP  arrays)  are  printed  by  EDITCR.  This  subroutine  is  called 
from  SETUP  and  at  the  end  of  every  run. 

EOSXPL :  EOSXPL  is  called  to  compute  pressures  and  the 

constant-energy  compressibilities  for  ideal  gasses  and 
explosives.  Gamma-law  equation-of -state  constants  are  given 
(Table.  2)  for  eighteen  explosives.  If  calculating  the 
pressure  for  cells  containing  explosive  material,  the  equa¬ 
tion  of  state  is  modified  for  cells  that  are  either  cut  by 
or  are  ahead  of  the  detonation  front.  The  modifications  to 
the  equation  of  state  are  such  as  to  give  zero  pressures 
ahead  of  the  detonation  front  without  breaking  down  the 
pressure  iteration.  See  Section  2.5  for  further  details. 
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An  entry  point  labeled  SETXPL  is  also  included  in 
EOSXPL.  SETXPL  is  called  in  setup  and  returns  the  gamma, 
detonation  velocity,  initial  density,  and  specific  available 
energy  for  the  explosive  being  set  up. 

EQST:  The  equation  of  state  constants  for  nineteen 

materials  are  stored  in  DATA  statements  in  EQST.  (These 
constants  are  displayed  in  Table  1.)  Given  a  code  number 
for  specifying  the  material,  a  density  and  an  energy,  EQST 
computes  a  pressure  and  a  cons tant-energy  compressibility. 

A  more  detailed  discussion  of  the  equation  of  state  is  given 
in  Section  2.3.  EQST  is  called  by  CDT  to  compute  pressures 
of  pure  cells  as  well  as  of  materials  in  mixed  cells.  EQST 
is  also  called  from  APLYPR  when  iterating  to  find  the  density 
of  material  in  a  cell  which  is  to  have  a  pressure  boundary 
condition  applied  to  it. 

ERROR:  This  subroutine  is  called  in  the  event  that 

certain  error  conditions  are  violated,  e.g.,  the  energy  sum 
error  exceeding  the  specified  limit.  ERROR  prints  a  message 
identifying  the  general  location  of  the  error  condition, 
lists  the  Z-block  variables  (first  150  words  of  blank  common), 
calls  EDIT  to  do  a  long  print  and  tape  dump,  and  then  calls 
exit. 


FA KM AT :  If  a  pressure  boundary  condition  is  being 
applied,  FAKMAT  identifies  all  interfaces  between  void  cells 
and  cells  which  have  had  their  pressures  defined  by  APLYPR. 
The  interface  pressures  are  redefined  in  a  manner  that  allows 
these  interfaces  to  be  treated  as  if  they  were  transmittive 
boundaries.  The  total  theoretical  energy  is  then  updated  to 
reflect  the  work  done  at  the  interfaces. 


88 


FLUX :  The  mass  and  volume  fluxes  across  the  right 

and  top  boundaries  of  a  mixed  cell  for  each  material  in  that 
cell  are  computed  in  FLUX.  These  fluxes  are  stored  in  the 
SAMPY,  SAMMP ,  VOLRT  and  VOLTP  arrays  and  executed  in  subroutine 
PH2  which  does  the  actual  transport  of  material  across  all 
cell  boundaries.  FLUX  uses  the  position  of  the  material 
interfaces  to  compute  the  fluxes  of  the  various  materials  in 
the  cell.  See  Section  4.5  for  a  more  complete  description  of 
transport  from  mixed  cells. 

INFACE :  The  tracer  particles  that  circumscribe  each 

material  package  are  moved  by  calling  MOVTCR  from  INFACE. 
Subroutine  CALFRC  is  then  called  by  INFACE  to  determine 
where  the  interfaces  cut  cell  boundaries  and  to  compute  the 
fractional  areas  for  each  material  at  these  interfaces.  For 
a  given  cycle,  INFACE  senses  which  cells  become  mixed  and 
which  become  pure.  INFACE  also  calls  FLUX  to  determine  the 
mass  and  volume  fluxes  across  mixed  cell  interfaces  and  calls 
ADJFLX  to  initialize  the  flux  arrays  and  to  set  fluxes  to 
remove  material  from  mixed  cells  which  become  pure. 

INPUT:  Instructions  for  running  problems  are  inter¬ 

preted  by  INPUT  which  can  either  start  or  restart  a  calculation. 
It  calls  SETUP  and  CARDS,  as  necessary,  to  prescribe  the  initial 
conditions  and  to  read  the  input  deck. 

LOCI J :  This  routine  is  used  to  identify  the  column  or 
row  location  of  a  tracer  point  whose  coordinates  are  in  centi¬ 
meter  units.  SETUP  calls  LOCIJ  when  converting  tracer 
coordinates  to  cell  units. 

MAP :  This  subroutine  is  called  by  EDIT  when  the  input 

parameter  MAPS  is  non- zero.  It  displays  the  properties  of 
each,  cell  in  the  active  grid  using  an  alphabetic  scale.  [On 
cycle  0  it  presents  the  entire  (IMAX  by  UMAX)  grid.]  One 
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obtains  contour  maps  of  the  compression,  pressure,  radial 
and  axial  velocities,  and  internal  energy.  Since  the  compression 
of  mixed  cells  is  not  computed,  asterisks  are  displayed  for 
mixed  cells  in  the  compression  map.  The  scale  of  each  map  is 
adjusted  according  to  the  current  minimum  and  maximum  values 
of  each  property. 

MOVTCR Subroutine  MOVTCR  moves  the  interface  and 
interior  tracers  by  looking  at  a  control  area  surrounding  each 
tracer  and  finding  which  cells  this  control  area  overlaps  and 
the  overlap  areas.  A  density-area  weighted  velocity  is  then 
computed  using  the  densities  and  velocities  of  the  cells 
overlapped  by  the  control  area.  Each  tracer  is  then  moved 
with  the  resulting  average  velocity. 

NEWMIX :  When  a  pure  cell  K  becomes  mixed,  NEWMIX 
assigns  to  it  a  location  M  in  the  mixed  cell  arrays  (SIE, 

XMASS,  RHO ,  etc.)  The  flag  of  cell  K  is  redefined  so  as  to 
associate  the  K  and  the  M  indices  (MFLAG(K)  =  M  +  100). 

The  fact  that  MFLAG(K)  is  greater  than  100  indicates  cell  K 
is.  a  mixed  cell.  Otherwise,  if  cell  K  were  pure,  MFLAG(K) 
would  be  equal  to  the  number  of  the  package  that  contained 
cell  K. 

NEWMIX  is  called  from  CALFRC  (called  from  INFACE  which 
is  usually  subcycled)  and  from  SETUP.  If  cell  K  becomes 
mixed  after  the  first  subcycle  of  INFACE,  NEWMIX  calls  FLUX 
to  bring  the  definition  of  the  flux  variables  for  cell  K 
to  the  current  time.  (See  Section  4.6  for  further  discussion 
of  mixed  cells . ) 

NEWR110 :  When  the  material  N  interface  first  enters 
a  cell,  CALFRC  calls  subroutine  NEWRUO  to  define  the  density 
of  material  N  for  cell  K  by  looking  at  the  density  of 
material  N  in  the  neighboring  cells  closest  to  the  interface. 
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(This  density  value,  stored  in  the  RIIO  array,  is  updated  in 
CUT  when  pressures  arc  equilibrated  at  the  beginning  of  each 
cycle.  If  the  cell  is  a  free  surface  cell,  the  densities  are 
updated  in  PHI  by  looking  at  the  velocity  field  and  in  Pl-12  by 
keeping  track  of  mass  and  volume  fluxes.) 

PHI :  The  effect  of  the  pressures  and  deviator  stresses 
in  updating  the  velocities  and  internal  energies  is  computed 
in  PHI.  [See  Section  2.2.2.]  The  effect  of  pressures  applied 
at  material  boundaries  and  the  updating  of  the  mass  densities 
in  free  surface  cells  are  also  handled  in  PHI. 

PH2 :  Mass  transport  and  the  associated  flux  of  momen¬ 
tum  and  energy  are  accounted  for  in  PH2 .  After  the  transport 
has  been  completed,  PH2  corrects  for  negative  mass  quantities 
in  mixed  cells  by  calling  subroutine  ZROMAS.  Before  printing 
the  symbolic  map  which  displays  the  material  package  numbers 
and  mixed  cell  locations,  PH2  redefines  the  flags  of  cells 
that  have  become  pure. 

PH3 :  The  cell-centered  deviator  stresses  are  updated 

each  cycle  in  PH3 .  (The  effect  of  these  deviator  stresses  on 
velocities  and  internal  energies  is  computed  in  PHI.)  PH3 
calls  STRNG  to  compute  the  yield  strength  of,  the  material  in 
a  cell.  If  the.  effects  of  strength  are  to  be  omitted 
(CYCPH3  £  0.),  PH3  insures  that  all  the  deviator  stresses 
are  zero  and  returns  to  PHI. 

PRESBC :  If  a  pressure  boundary  condition  is  desired, 

the  user  must  supply  this  subroutine  which  computes  a  pressure, 
density  and  energy  for  any  given  time  and  location  in  the  grid. 
See  Section  6  for  a  more  detailed  explanation. 

PRSV1S:  Subroutine  PRSVIS  is  called  from  Pill  and  com¬ 
putes  a  pseudo- viscous  pressure  which  is  added  to  the  interface 
pressure  terms.  This  artificial  viscosity  helps  minimize  in¬ 
stabilities  in  low  velocity  stagnation  regions. 


RMVMAT :  If  a  pressure  boundary  condition  is  being 

used,  RMVMAT  (which  is  called  whenever  a  problem  is  being 
restarted  and  at  the  end  of  PHI)  removes  any  mass  and  energy 
which  was  added  by  APLYPR  and  updates  the  total  theoretical 
energy. 

SETUP :  This  routine  generates  the  shaped  charge  prob¬ 

lems  to  be  run.  It  defines  the  X,  Y,  DX,  DY  and  TAU  arrays 
for  variable  as  well  as  constant  zoning  in  either  cylindrical 
or  rectangular  geometries.  The  properties  of  material  packages 
are  determined  by  reading  the  INPUT  cards.  SETUP1  is  called 
to  fill  the  cells  with  the  necessary  material,  SETUP4  is  called 
to  place  the  interface  tracers,  and  SETUP5  is  called  to  place 
the  interior  tracers.  The  total  theoretical  energy  is  set 
equal  to  the  total  energy  in  the  problem  and  a  tape  dump  is 
written. 

SETUP  1 :  After  SETUP  has  read  the  INPUT  cards,  SETUP1 
is  called.  SETUP1  first  checks  to  make  sure  that  the  INPUT 
data  defines  an  acceptable  shaped  charge  configuration  and 
then  calls  SETUP2  to  determine  MFLAG  and  compute  the  mass  in 
every  cell.  SETUP1.  then  sets  the  U,  V,  and  AIX  arrays  for 
each  cell. 

SETUP2:  SETUP 2  sets  the  MFLAG,  AMX ,  XMASS,  and  RHO 
arrays  for  each  cell  in  the  grid.  The  mass  for  each  material 
is  computed  by  calling  SETUP3  to  determine  the  partial  volume 
each  material  occupies  in  a  cell. 

SETUPS :  SETUP3  determines  if  a  cell  is  void,  pure,  or 

mixed  and  returns  the  partial  volume  occupied  by  each  material 
and  the  void  in  the  cell. 

SETUP 4 :  The  interface  tracers  (TX,  and  TY  arrays)  are 
placed  around  the  material  packages  and  the  void  by  SETUP 4. 
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SETUP 5 :  The  passive  interior  tracers  (XP  and  YP  arrays) 
are  placed  by  subroutine  SETUP5. 

STRNG :  The  yield  strength  of  the  material  in  a  cell 

is  computed  by  STRNG.  (The  strength  constants  for  the  liner 
material  are  defined  by  input  cards  when  the  problem  is 
generated.)  The  yield  strength  of  the  material  in  a  mixed 
cell  is  a  volume  weighted  average  of  the  yield  strengths  of 
the  various  materials  in  the  cell*  STRNG  is  called  from  PH5. 

Z ROMAS :  If  PH2  transports  more  mass  out  of  a  mixed 
cell  than  is  in  the  cell,  ZROMAS  is  called.  ZROMAS  determines 
which  cells  received  the  excess  mass.  Enough  material  is 
removed  from  these  cells  to  bring  the  total  mass,  momentum 
and  energy  of  the  over-evacuated  material  in  the  cell  back  to 
zero. 


93 


8.3  GENERATING  PROBLEMS 


8.3.1  General  Capabilities  and  Limitations 

The  procedure  for  setting  up  a  problem  to  be  run  by 
BRLSC  has  been  simplified  to  allow  problems  to  be  defined  as 
easily  as  possible.  Many  of  the  variables  have  predetermined 
default  values  assigned  to  them  so  that  it  is  not  necessary 
to  enter  them  unless  they  are  to  be  changed. 

8 . 3 . 1 . 1  Modes 

BRLSC  is  capable  of  running  problems  in  any  of  three 
basic  modes.  They  are: 

1.  In  Mode  1,  complete  shaped  charge  problems 
are  set  up  and  run.  Both  the  explosive  and 
the  liner  are  in  the  grid  and  the  problem 
can  be  run  as  far  in  time  as  desired  with  no 
additional  input.  See  Section  7.1  for  an 
example  of  a  problem  run  in  Mode  1. 

2.  In  Mode  2,  only  the  liner  is  set  up  in  the 
grid.  The  explosive  driving  pressures  are 
applied  along  the  outside  surface  of  the 
liner  in  the  form  of  a  time  and  position 
dependent  pressure  boundary  condition.  The 
pressure,  density  and  specific  internal 
energy  of  the  driving  explosive  are  entered 
into  the  calculation  each  cycle  from  subroutine 
PRESBC.  The  user  must  supply  this  subroutine 
for  Mode  2  calcul ati ons .  The  values  returned 
can  be  determined  in  any  way  the  user  wishes: 
analytical  equations,  table  lookup,  another 
calculation  (Eulerian  or  Lagrangian) ,  etc. 

The  resulting  liner  motion  is  then  calculated. 

See  Section  7.2  for  an  example  of  Mode  2. 
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3.  In  Mode  5,  the  liner  only  is  set  up  and  given 
an  initial  velocity.  Mode  3  can  be  used  to 
investigate  such  problems  as  an  initially 
free  flying  plate  being  driven  into  the  plane  ' 
of  symmetry.  See  Section  7.3  for  an  example 
of  Mode  3. 

8. 3. 1.2  Restrictions  and  Limitations 

When  setting  up . a  problem,  several  restrictions  and 
limitations  must  be  kept  in  mind.  These  are: 

1.  No  more  than  two  material  packages  may  be  set 
up,  one  for  the  liner  and,  if  running  in  Mode 
1,  one  for  the  explosive. 

2.  All  the  material  packages  must  be  entirely 
contained  in  the  grid. 

3.  The  half  angle  of  the  liner  must  be  between 
five  and  eighty-five  degrees. 

4.  If  the  liner  is  set  up  with  a  round  tip,-  the 
inside  radius  of  the  tip  must  be  less  than  the 
inside  radius  of  the  base  of  the  liner. 

5.  When  setting  up  a  problem  in  Modes  1  or  2 ,  the 
initial  radial  velocity  must  be  zero.. 

6.  The  left-hand  boundary  of  the  grid  (i.e.,  the 
axis)  is  always  considered  to  be  reflective. 

The  other  three  boundaries  are  always  trans- 
mittive.  . 

7.  In. Mode  1,  the  top  of  the  explosive  package 
must  be  above  the  tip  of  the  liner,  the 
radius  of  the  explosive  package  must  be 
greater  or  equal  to  the  outside  radius  of 
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the  base  of  the  liner,  and  the  axial  position 
of  the  bottom  of  the  explosive  is  forced  to 
be  equal  to  that  of  the  liner. 

8.  The  diagonal  dimension  of  every  cell  which 
contains  liner  material  should  be  less  than 
the  thickness  of  the  liner.  Best  results  are 
obtained  if  all  cells  initially  containing  liner 
material  are  square,  i.e.,  DX  =  DY. 

9.  It  is  suggested  that  when  setting  up  a  problem 
in  Modes  1  or  2,  that  the  initial  axial  velo¬ 
city  be  zero. 

8.3. 1.3  Capabilities 

The  capabilities  of  BRLSC  for  setting  up  and  running 
shaped  charge  problems  are  generalized  below. 

1.  Choice  of  operating  mode: 

•  Two  materials,  i.e.,  liner  and  explosive 

•  Liner  and  pressure  boundary  condition 

•  Liner  only 

2.  Choice  of  liner  material: 

•  Tungsten 

•  Copper 

•  Iron 

•  Aluminum 

•  Beryllium 

•  Titanium 

•  Nickle 

•  Molybdenum 

•  Thorium 

•  Lead 
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3.  Choice  of  Explosives: 


o  Baratol 
9  Comp  B ,  Grade  A 
o  Comp  B-3 

•  Cyclotol  (77/23) 
o  Datb 

o  HMX 
o  LX- 01 -  0 
o  LX- 04-1 
©  LX-07-0 

9  Nitroglycerine 

•  Nitromethane 

e  Octol  (77.6/22.4) 

•  PBX-90 10 

•  PBX-9011 

9  PBX-9404-03 

•  PETW 

•  RDX 

o  TNT  • 

4.  Choice  of. liner  half  angle,  thickness,  base 
radius,  and  tip  geometry  (round  or  pointed). 

5.  Choice  of  explosive  dimensions. 

6.  Choice  of  spherical  or  plane  detonation  front 
If  spherical,  choice  of  initial  radius  of 
curvature . 

7.  Option  of  including  or  not  including  the 
effects  of  material  strength  in  the 
calculation. 
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8.3.2 


Format  of  Input  Cards 


A  BRLSC  problem  is  generated  by  three  kinds  of  input 

data: 

1.  The  identification  header  card  read  by  INPUT 
that  contains  alphanumeric  symbols  in 
columns  2  through  72. 

2.  The  input  data  read  by  subroutine  CARDS 
that  defines  variables  in  the  Z-block  (the 
first  ISO  words  of  blank  common) . 

3.  The  input  data  read  by  subroutine  SETUP 
that  defines  the  DX ,  DY  arrays  when 
either  or  both  are  variable  and  the 
boundaries  and  properties  of  the  material 
packages. 

The  input  cards  read  by  CARDS  use  the  following  format: 

(II)  (E9.4)  (E9.4) 

Col.  7  Cols.  8-16  Cols.  17-35  etc. 

N  Z(M)  Z (M+l) 

The  value  of  L  for  each  input  variable  is  listed  in  the 
discussion  that  follows.  When  L=0  the  variable  defined  is 
typed  real.  When  L=2,  the  variable  defined  is  typed  integer. 
When  L=1  the  variable  defined  is  typed  real  and  the  card 
is  the  last  one  to  be  read  until  CARDS  is  called  again.  The 

t 

value  of  M  is  the  location  in  blank  common  of  the  variable 
being  defined  in  Cols.  8-16.;  this  location  number  is  also 
listed  below  for  each  input  variable.  The  value  of  N  indi¬ 
cates  how  many  consecutively  stored  variables  of  the  same  • 
type  are  defined  by  a  card  (usually  N=l) .  Z(M)  must  be 
punched  with  an  E  format  even  when  it  is  defining  an  inte¬ 
ger  variable. 


(ID  (15) 

Col.  1  Cols  .  2-6 
L  M 


98 


The  format  of  the  cards  read  by  SETUP  varies  and  is 
described  in  the  following  discussion  of  the  input  variables. 

8.3.3  Description  of  INPUT  Variables 

The  following  list  of  INPUT  variables  is  ordered  as 
the  BRLSC  setup  INPUT  deck  should  be  ordered.  These  variables 
are  broken  into  four  groups.  The  first  three  groups  of  cards 
are  read  by  CARDS.  The  fourth  group  .is  read  by  SETUP  and  the 
last  dummy  card  is’  read  by  CARDS. 

The  first  group  contains  the  cards  which  must  be  de¬ 
fined  in  order  to  run  a  problem.  The  second  group  lists  the 
variables  which  have  default  values  assigned  them  and  which 
are  generally  omitted  unless  these  default  values  are  to  be 
changed.  The  third  group  contains  variables  which  are  either 
calculated  or  zero  at  setup  time.  Except  in  rare  instances 
when  restarting  a  problem  (i.e.,  switching  to  a  new  dump 
tape,  etc.),  there  is  no  need  to  ever  input  or  change  them. 

8 . 3 . 3 . 1  Group  1  (Most  Common  Variables  Read  by  CARDS) 


Variable 

.  Column 

Location  in 

Name 

One  Flag 

Blank  Common 

Definition 

PK(1) 

1 

151 

Problem  number.  Any  number 
from  .0001  to  99.9999. 

Should  be  identical  to  PROB. 
This  is  the  second  input  card 
of  a  restart  as  well  as  a 
setup  deck.  (The  first  card 
is  the  identification  header 
card.) 

PROB 

1 

Problem  number  [same  as 
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Variable  Column  Location  in 
Name  One  Flag  Blank  Common 

NFRELP  2  5 


NDUMP7  2  6 

ICSTOP  2  7 


IGM  2  21 

IMAX  2  33 

JMAX  2  35 

MAPS  2  42 


PRDHLT  45 


Definition 

Gives  frequency  of  "long" 

EDIT  prints  which  gives  velo¬ 
cities,  energy,  compression 
and  stresses  for  all  cells 
in  active  grid.  A  "short" 

EDIT  gives  this  information 
only  for  the  axis  column. 

When  NFRELP=2,  every  second 
EDIT  print  is  "long". 

Gives  frequency  of  restart 
tape  dumps.  Program  will 
dump  only  when  EDIT  prints. 

If  NDUMP7  =  5,  a  restart  .dump 
will  be  made  every  5th  time 
EDIT  prints. 

Gives  cycle  at  which  execu¬ 
tion  will  stop  if  the  user 
wants  to  specify  a  cycle 
rather  than  a  time  on  which 
to  stop.  If  stopping  on 
time  omit  this  card. 

When  IGM=0.  code  uses  cylin¬ 
drical  coordinates.  When 
IGM=1.,  code  uses  plane 
coordinates . 

The  number  of  columns  in  the 
grid. 

The  number  of  rows  in  the 
grid. 

When  MAPS  >  0.,  part  of  the 
EDIT  print  is  a  set  of  sym¬ 
bolic  maps  of  the  compression, 
pressures,  velocities,  and 
specified  internal  energies 
of  cells  in  the  active  grid. 

The  time  (sec)  between  EDIT 
prints.  If  printing  on 
cycle  intervals  omit  this 
card . 
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Variable  Column  Location  in 

Name  One  Flag  Blank  Common  Definition 


IPCYCL  2 


CYCPH3 


NMXCLS  2 


NTPMX  2 


INTER  2 


49 


70 


73 


78 
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The  number  of  cycles  between 
EDIT  prints.  If  printing  on 
time  intervals  omit  this  card. 

The  number  of  passes  through 
subroutine  PH3  (strength 
phase)  in  one  time  step. 

Usually  CYCPH5  =  2.,  but 
larger  values  can  help  to 
dampen  os cillati ons  .  Omit 
this  card  when  the  calculation 
is  pure  hydro  (no  strength 
effects)  . 

The  maximum  number  of  mixed 
cells  to  exist  in  the  grid  at 
any  one  time  during  the  cal¬ 
culation.  Consult  the  dimen¬ 
sions  of  the  mixed  cell  arrays 
in  MX CELL  COMMON  and  BLANK 
COMMON. 

Used  in  SETUP  to  compute  the 
number  of  tracers  to  be  placed 
around  each  material  package. 
The  code  will  place  the  tracers 
in  such  a  way  that  no  package 
will  contain  more  than  NTPMX 
tracers.  This  maximum  cannot 
be  greater  than  the  dimensions 
of  the  TX  and  TY  arrays. 

A  special  editing  flag  for 
debugging.  When  INTER  =1., 
EDIT  prints  after  PHI  as  well 
as  after  CDT  on  print  cycles. 
When  INTER  =7.,  PI  1 2  prints 
energy  totals  after  updating 
each  cell.  The  later  option 
is  in  addition  to  the  extra 
EDIT  prints  and  involve  many 
pages  of  output. 


Variable  Column  Location  in 
Name  One  Flag  Blank  Common 


Definition 


NODUMP  2 


DXF 


DYF 


TSTOP  1 


96  When  NODUMP  =  1.,  EDIT  will 

not  write  any  tape  dumps. 
(This  flag  overrides  the 
NDUMP7  option,  but  does  not 
prevent  SETUP  from  writing 
the  cycle  0  dump.) 

136  The  x-dimension  of  cells  (cm) 
when  the  grid  has  constant 
zoning  in  the  x-direction. 
Omit  this  card  if  the  x- 
dimension  of  cells  varies. 

137  The  y-dimension  of  cells 
(cm)  when  the  grid  has 
cons tant  zoning  in  the  y- 
direction .  Omit  this  card 
if  the  y-dimension  of  cells 
varies. 

SO  The  time  (sec)  at  which  the 

calculation  will  stop.  This 
card  (with  a  "1"  in  column 
one)  should  not  be  omitted 
even  when  stopping  on  a 
specified  cycle,  and  should 
always  be  the  last  card  of 
a  restart  deck. 
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8.3. 3. 2  Group  2  (Variables  with  Default  Values  Assigned) 


Variable  Column  Location  in 

Name  One  Flag  Blank  Common  Definition 


IPR 


PRCNT 


DM  IN 


15 


16 


24 


Maximum  number  of  iterations 
CDT  will  perform  when  attemp¬ 
ting  to  equilibrate  the 
pressures  of  materials  in  a 
mixed  cell.  If  the  pressures 
are  not  within  an  epsilon 
(PRCNT)  of  the  average 
pressure  after  IPR  iterations, 
an  error  exit  occurs. 

(DEFAULT  =  30) 

Convergence  minimum  for  the 
pressure  equilibration  of 
materials  in  a  mixed  cell. 

All  pressures,  Pj_,  must 
satisfy  the  following: 

l(Pav-pi)/Pavl<PRCNT- 
(DEFAULT= . 001) 

The  maximum  relative  -error 
in  the  energy  sum  that  the 
user  wants  to  tolerate. 

The  error  is  computed  as 
follows:  . 

KMAX 

E 

K=  1 

where  E,  is  total  energy 
of  cell  KK  and  ETH  is 
theoretical  energy  in  the 
gri d- -  computed  in  SETUP  and 
updated  in  PHI,  PI 1 2  ,  APLYPR 
and  RMVMAT .  (DEFAULT  =  10“ 3) 


Ek  -  ETH 


/ETH 


/ 
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Variable  Column  Location  in 

Name  One  Flag  Blank  Common  Definition 


CYCMX  69  The  number  of  passes  through 

subroutine  INFACE  in  one  time 
step  to  minimize  transport 
noise  near  interfaces. 


Usually  CYCMX  =  2 ,  but  when 
calculation  costs  permit, 

CYCMX  =  8  is  a  good  maximum. 
(DEFAULT  =  2) 

PM  IN 

86 

A  pressure  cutoff  (DEFAULT 
=  5  x  106)  . 

ROEPS 

110 

A  round-off  epsilon  used  to 
define  cutoffs  in  the  cal¬ 
culation.  (DEFAULT  =  10" 5) 

FINAL 

113 

The  final  value  of  the 
stability  fraction  used 
in  determining  the  time 
step.  (See  STAB) (DEFAULT  =  .4) 

STAB  • 

139 

The  initial  value  of  the 
"stability  fraction"  for 
the  calculation  of  DT.  If 
FINAL=0.,  STAB  is  constant. 
Otherwise  its  value  increases 
from  the  input  value  to 

FINAL  in  a  geometric  pro¬ 
gression.  (DEFAULT  =  .001) 

DTMIN 

144 

The  minimum  value  for. the 
time  step.  The  program 
gives  an  error  exit  if  CDT 
calculates  a  At  <  DTMIN. 
(DEFAULT  =  10 -1  *) 

BBAR 

149 

A  constant  used  in  calcula- 

tion  of  local  sound  speed  of 
all  materials  except  ideal__ 
gases.  C  =  C  +  BBAR  •  (/P) . 
(DEFAULT  =  .5j 
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8.3. 3.3  Group  3  (Seldom  Used  Variables) 


Variable  Column  Location  in 

Name  One  Flag  Blank  Common  Definition 


KUNITR  2  14 


KUNITW  2  17 


CVIS  27 


NUMSCA  2  43 


PRLIM  44 


The  tape  unit  INPUT  will  read 
from.  KUNITR  will  be  set 
equal  to  7  unless  defined  by 
this  input  card. 

The  tape  unit  EDIT  and  SETUP 
will  write  on.  KUNITW  will 
be  set  equal  to  7  unless 
defined  by  this  input  card. 

A  .flag  that  describes  the 
boundary  condition  at  the 
bottom  of  the  grid.  If 
CVIS=0.,  the  boundary  will 
be  reflective.  If  CVIS=-1., 
the  boundary  will  be 
transmittive . 

The  number  of  times  the  cycle 
or  time  interval  between  EDIT 
prints  is  increased.  (See 
PRLIM) 

The  time  or  cycle  at  which 
the  EDIT  print  frequency 
is  changed.  Example:  you 
wish  to  print  every  10“ 8  sec 
until  T  =  10" 7  sec;  and  every 
10-7  sec  until  T  =  10~6  sec', 
and  every  10"6  sec  thereafter. 

Set:  PRDELT  =  10" 8 

IPCYCL  =0. 

PRLIM  =  10" 7 

PRFACT  =10 

NUMSCA  =2. 

NOTE:  When  PRDELT  is  increased 
by  a  factor,  PRLIM  is  also  in¬ 
creased  by  t lie  same  factor 
for  the  next  rescaling.  If 
you  want  a  constant  print 
interval  omit  NUMSCA,  PRLIM, 
PRFACT. 


105 


Variable  Column  Location  in 

Name  One  Flag  Blank  Common  Definition 


PRFACT 

46 

The  factor  by  which  the  print 
interval  (PRDELT  or  IPCYCL) 
and  tlie  print  limit  (PRLIM) 
are  increased.  (See  PRLIM) 

11 

2 

47 

The  number  of  columns  in  the 
grid  that  have  non-zero  velo¬ 
cities  or  energies  plus  2. 

11  and  T2  define  the  "active" 
grid.  (Most  calculations  are 
done  inside  the  active  grid.) 

12 

2 

48 

The  number  of  rows  in  the  grid 
that  have  non-zero  velocities 
or  energy  plus  2 . 

IVARDY 

2 

54 

When  the  grid  is  variably 
zoned  in  the  y-direction 
(DX’s  not  constant), 

IVARDY* 1.  When  the  grid  has 
constant  zoning  in  the  y- 
direction  (DYF>0.),  IVARDY* 0 . 

N6 

2 

56 

This  parameter  causes  the 
program  to  set  negative 
pressures  to  zero  in  the 
rows  of  the  grid  below 

J  =  N6.  Omit  this  card  if 
you  want  to  allow  negative 
pressures  everywhere  in  the 
grid.  If  you  want  to  set 
negative  pressures  to  zero 
everywhere,  set  N6  =  104. 

GAMMA 

62 

(a+1.)  in  gainma-law  equation 
of  state.  Used  only  for  an 
ideal  gas. 

NMAT 

2 

6  8  , 

The  number  of  material 
packages,  exluding  the  void 
package.  The  maximum  number 
of  packages  is  two. 
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Variable 

Name 


Column  Location  in 
One  Flag  Blank  Common 


Definition 


IVARDX  2  83 


EMIN  85 


DTFIX  143 


When  the  grid  is  variably 
zoned  in  the  x-dircction  (DX 1 s 
not  constant),  IVARDX=1. 

When  the  grid  has  constant 
zoning  in  the  x-direction 
(DXF>0),  IVARDX- 0 . 

The  minimum  specific  internal 
energy  of  an  ideal  gas  to  be 
used  in  the  pressure  iteration. 
Usually  EMIN  =  106. 

Gives  the  maximum  DT  to  be 
used  in  the  calculation. 

If  DTFIX>DMIN ,  then 
DT=min(DT, DTFIX)  where  DT 
is  computed  each  cycle  by 
CDT .  Omit  this  card  if  no 
maximum  is  to  be  place  on 
the  time  step. 


8. 3. 3. 4  Group  4  (Input  Read  by  SETUP) 

1.  Defining  Coll  Dimensions  IVhcn  Variable,  User  can 
specify  variable  DY  and/or  DX  values.  In  the  case  that  both 
are  variable,  the  DY  array  is  defined  first.  Both  arrays  are 
read  using  the  following  format: 


Variable  Name 

Columns 

Definition 

NT  (I)  ,  1  =  1,4) 

1-4 

NT(I)  is  the  number  of  zones  that 

5-8 

have  dimension  TEMP (I).  When  all 

9-12 

DY T s  are  defined,  set  next  NT-999. 

13-16 

Likewise  for  DX's. 

TEMP (I) ,1=1,4 

21-30 

NOTE:  When  DY’s  are  variable, 

31-40 

exactly  JMAX  DYfs  must  be  defined. 

41-50 

Likewise,  exactly  I MAX  DX's  must  be 

51-60 

defined . 

2. 

Def inin 

g  the  Liner. 

i.  —  ■ 

Variable 

Card 

Columns 

Name 

Number 

(Format) 

Definition 

MODE 

1 

1 

(ID 

If  MODE-1,  both  the  liner  and  ex¬ 
plosive  are  set  up.  If  M0DE=2, 

only  the  liner  is  set  up  and  a 
pressure  boundary  condition  will 
be  applied  along  the  outer  free 
surface.  If  M0DE=5,  only  the  liner 
will  be  set  up  with  some  initial 
velocity.  See  Section  S.3.1.1  for 
details . 


MAT(l)  2  1-5  Material  code  number  for  the  liner 

(15)  1  <  MAT(l)  <  10.  (i.e.,  if  MAT(1)=2, 

the  liner  is  copper.)  Sec  Table  1, 
Section  2.3.1  for  list  of  liner 
materials  available. 
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Variable 

Name 

Card 

Number 

Columns 

(Format) 

Definition 

NTHCK 

2 

6-10 

(15) 

Indicates  the  number  of  interior 
tracers  to  be  placed  across  the 
thickness  of  the  liner.  If  no 
interior  tracers  are  desired, 
NTHCK=0 . 

ALPHA 

2 

11-20 
(E10  .2) 

ALPHA  is  the  half  angle  of  the  top 
of  the  liner.  5°  <_  ALPHA  <_  85°. 

YTIPI 

2 

21-30 
(E10. 2) 

YTIPI  is  the  Y  coordinate  of  the 
inside  tip  of  the  liner.  If 

YTIPI  >  0.,  the  tip  is  placed 
YTIPI  cm  from  bottom  of  grid. 

If  YTIPI  <  0,  the  tip  is  placed 
at  | YTIPI |  cells  from  the  bottom 


of  the  grid. 

RTIPI 

2 

31-40 
(E10 .2) 

Radius  in  cm  of  the  inside  of  the 
tip  of  the  liner.  If  RTIPI  <  0, 
the  liner  is  set  up  with  a  pointed 
tip  . 

THICK 

2 

41-50 

(E10.2) 

Thickness  (cm)  of  the  liner. 

RLINI 

2 

51-60 
(E10. 2) 

Radius  (cm)  of  the  inside  of  the 
base  of  the  liner. 

UUR(l) 

3 

1-10 
(E10. 2) 

Initial  radial  velocity  of  liner 
(cm/sec).  UUR(l)  will  be  set 
=0.  in  Modes  1  and  2. 

WA  ( 1 ) 

5 

11-20 

(E10.2) 

Initial  axial  velocity  (cm/sec) 
of  entire  shaped  charge  (usually 
=0.  in  Modes  1  or  2). 

SSIEN(l) 

3 

21-30 

(E10.2) 

Initial  specific  energy  (ergs/g)  of 
liner  material. 

RHOIN(l) 

3 

31-40 
- (E10 .2) 

Initial  density  (g/cm3)  of  liner 
material . 
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Variable 

Card 

Columns 

Name 

Number 

(Format) 

Definition 

CZERO(l) 

4 

1-10 

\ 

Y 

|  used  to  determine  yield 

(E10.2) 

0 

|  strength  of  the 

liner  material . 

STK1(1) 

4 

11-20 

C 

K  Y.S.  =  (y  +  C*P 

)(  1  -E/E  ) where 

(E10.2) 

P  =  pressure  and 
|  internal  energy. 

L  =  specific 

STEZ(l) 

4 

21-30 

E 

(E10.2) 

0 

J 

> 

RMU(l) 

4 

31-40 

Rigidity  modulus  for 

liner  material . 

(E10.2) 

AM DM (1) 

4 

41-50 

Failure  criterion  for 

liner  material 

(E10.2) 

If 

the  compression  of 

the  liner 

material  in  a  cell  is  less  than 
AM DM ,  the  cell  has  failed  and  all 
tensile  stresses  are  set  to  zero. 


The  following  card  is  used  only  if. MODE  =  1.  It  describes  the 
explosive  package. 


M ATX PL 

5 

1-5 

(IS) 

The  absolute  value  of  MATXPL  de - 
fines  tlie  type  of  explosive  to  be 
used.  (See  Table  2,  Section  2.3.2 
for  list  of  available  explosives) .  If 
MATXPL  is  >  0 ,  a  plane  detonation 
front  will  be  used.  If  MATXPL  <  0, 

a  curved  detonation  front  is  used. 

XPLRAD 

5 

11-20 

(E10.2) 

If  the  problem  is  being  set  up  to 
run  with  a  plane  detonation  front, 
XPLRAD  will  be  set  =0.  If  running 
with  a  curved  detonation  front, 

XPLRAD  is  the  initial  radius  (cm) 
of  the  detonation  front.  (If 

XPLRAD  =  5.,  the  explosive  was 
detonated  at  a  point  5  cm  above 
the  top  of  the  explosive  package.) 

RITXPL 

5 

21-30 

(E10.2) 

The  radius  (cm)  of  the  outside  of 
the  explosive. 

no 


Variable  Card  Columns 

Name  Number  (Format)  Definition 


TOPXPL 

5 

51-40 

(E10.2) 

The  distance  (cm)  from  the  outside 
tip  of  the  liner  to  the  top  of  .the 
explosive . 

EMIN. 

5 

41-50 
(E10 .2) 

The  minimum  specific  internal  energy 
(ergs/g)  of  the  explosive  to  be  used 
in  the  pressure  iteration.  If  EMIN 
is  entered  as  0,  the  program  will 
give  it  a  value  of  10  . 

The  last  card  on  the  setup  deck  (read. by  CARDS)  is 


Variable  Column  Location  in 

Name  One  Flag  Blank  Common  Definition 


EMOB  1  150  Dummy  end  card  of  SETUP  deck. 

Omit  from  restart  deck. 

EMOB  =  0. 


Ill 


8.3.4  Sample  Input  Deck  for  Setting  up  a  BRLSC  Problem 


The  input  deck  for  setting  up  a  mode  1  shaped  charge 
problem  is  listed  below.  This  problem  has  been  set  up  and 
run  to  a  time  of  one  microsecond  on  the  version  of  BRLSC 
listed  in  this  report. 
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The  following  pages  show  part  of  the  printed  output 
from  the  sample  problem  whose  input  deck  is  described  above. 

The  calculation  was  for  a  shaped  charge  with  a  .1414  cm  thick 
copper  liner  and  a  Comp  B  explosive.  The  half-angle  of  the 
liner  was  21  degrees,  the  radius  of  the  base  of  the  liner 
was  .4  cm,  and  the  radius  of  the  tip  of  the  liner  was  .1414  cm. 
A  curved  detonation  front  was  employed  and  strength  effects 
were  ignored. 
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8.4  RESTARTING  PROBLEMS 


During  the  calculation,  BRLSC  periodically  writes  on  a 
tape  or  drum  file  the  problem  parameters  and  the  current  state 
of  the  material  in  each  cell.  By  reading  this  file  and  a  few 
input  cards,  the  code  can  "restart”  and  continue  a  calculation 
from  an  intermediate  point. 

The  three  cards  that  must  be  included  in  a  restart  deck 
are  as  follows : 

1.  the  heading  card 

2.  the  restart  card 

3.  the  last  card 

The  heading  card  can  contain  any  alphanumeric  symbols 
between  columns  2  and  72.  ■  The  "restart"  card  defines  the  four 
variables  described  below. 


Variable 

Name  Columns  Definition 


PK(1) 

8-16 

Problem  number. 

PK(2) 

17-25 

Restart  cycle  number. 

PK(3) 

26-34 

Restart  flag. 

=-l.  prints  long  edit  of  restart 
--2.  prints  short  edit  of  restart 

cycle 
cycle . 

PK  (4) 

35-43 

Number  of  material  packages  in  the  problem 
(excluding  the  void  package). 

This  restart  card  is  read  by  CARDS  and  columns  1  through  7 
are  punched  as  follows 


(Columns ) 


which  indicate  there  are  four  words  defined  by  the  card,  the 
first  one  stored  in  the  151st  location  of  blank  common. 

The  last  card  defines  TSTOP,  the  time  to  stop  the  cal¬ 
culation.  This  must  be  the  last  card  of  the  restart  cycle. 
The  TSTOP  card  is  punched  as  described  above  for  a  setup 
deck ,  e  .g  .  , 


1 

2  3  4  5  6 

7 

8  9  10  12  13  14  15  16 

(Columns) 

1 

.5  0  1 

1 

1.0  0  6 

Many  problem  parameters  stored  in  blank  common  can  be 
redefined  when  restarting  at  an  intermediate  point  in  the 
calculation.  For  example,  the  number  of  times  PH3  is  sub- 
cycled,  the  frequency  of  printing,  or  the  round-off  epsilon 
can  be  changed  by  including  in  the  restart  deck  cards  that 
redefine  these  variables.  These  extra  cards  must  be  added 
between  the  '’restart"  card  and  the  last  card,  i.e.,  between 
the  two  cards  that  have  a  "1"  in  column  1.  All  cards  except 
the  heading  card  of  the  restart  deck  are  read  by  the  CARDS 
subroutine  with  the  format  described  in  Section  8.3.2. 

Examples  of  Restart  Decks 

The  three  cards  listed  below  could  be  used  to  restart 
a  problem  from  cycle  100  and  run  until  5  usee.  There  are  two 
materials  in  the  problem  (Mode  1)  and  a  short  print  is  de¬ 
sired  on  the  restart  cycle. 

SHAPED  CHARGE  TEST  PROBLEM  11.11 


151 

4 

11.11 

100.0 

-2.0 

50 

T 

5.0 

-06 

1 

1 


2.0 


The  cards  listed  below  could  be  used  to  restart  a 
problem  from  Cycle  100,  and  run  to  Cycle  200.  There  is  only 
one  material  in  the  problem  (Mode  2  or  3)  and  a  long  print 
is  desired  on  the  restart  cycle.  The  print  material  will 
also  be  changed  so  that  a  long  edit  print,  with  symbolic  maps, 
will  be  made  every  25  cycles. 

SHAPED  CHARGE  TEST  PROBLEM  22.22 


1 

151 

4 

22.22 

100.0 

-1.0 

2 

5 

1 

1.0 

2 

7 

1 

200. 

2 

42 

1 

1.0 

2 

49 

1 

25.0 

1. 

50 

1 

0. 

/ 
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8.5  DICTIONARY  OF  IMPORTANT  VARIABLES 


This  section  includes  a  description  of  the  use,  units, 
dimensions,  and  locations  of  all  important  variables  in  named 
or  blank  commons  and  some  variables  local  to  the  subroutines. 

Location 

(NAME)  Variable  is  local  to  subroutine  NAME. 


NAME 

Variable  is  in  common  block  NAME. 

B.C. 

The  variable  is  in  Blank  Common  or  equi 
valenced  to  a  variable  in  Blank  Common. 

- 

=  Z(N) 

The  variable  is  equivalenced  to  Z  (N)  , 
the  first  array  in  Blank  Common.  These 
variables  are  usually  used  in  setting 
up  and  restarting  problems. 

Dimension 

IDX 

Number  of  cells  in  radial  direction 

(=40) 

JDX 

Number  of  cells  in  axial  direction 

(=120) 

NMX 

Number  of  material  packages 

(=2) 

MDX 

Number  of  mixed  cells 

(=400) 

NTPMX 

Number  of  tracers  per  package 

(=300) 

IDXX 

IDX+1 

(  =  41). 

JDXX 

JDX+1 

(=121) 

NMXX 

NMX+1 

(=3) 

JDX2 

2 -JDX 

(=240) 

KDX 

IDX • JDX+2 

(=4802) 

JDXP 

10 ‘NMX+JDX2 

(=260) 

MDX  MX 

NMXX • MDX +1 

(=1201) 

KDX4 

(IDX/2+1) • (JDX/ 2+1) 

(=1281) 

KDXP 

MAX [KDX, NMX* (JDX+5) ,  2*NMXX*MDX] 

(=4802) 
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Variable 

Name 

Location 

Dinens ion 

Units 

Definition 

AIX 

B.C. 

(KDX)  • 

ergs/g 

Specific  internal  energy  in  a  cell.  ( I MAX 
by  JMAX  array.) 

ALE 

(MAP) 

(41) 

*  " 

This  array  has  alphabetic  characteris  for 
positive  values  in  the  density,  velocity  and 
energy  maps.  Defined  in  a  DATA  statement. 

AM  DM 

B.C. 

(NMX) 

• 

INPUT  parameter.  A  material  i  with  com¬ 
pression  >  AMDM^  is  considered  solid. 

Usual  value:  0.9S  to  0.99. 

Used  in  CDT  in  testing  whether  to  allow 
negative  pressures  (tensions)  in  pure  cells. 
Used  in  PH3  to  calculated  SOLID  =  AM  DM  *  RHC2 
which  determines  whether  material  stresses 
are  nonzero. 

AMMP 

(PII2) 

g 

Mass  transported  across  right  boundary  of 
a  cell.  (Sec  Appendix  B) 

AMMU 

(PH2) 

-  - 

cm-g/sec 

Radial  momentum  transported  across  the  bottom 
boundary  of  a  cell.  (See  Appendix  B) 

AMMV 

(PH2) 

cm-g/scc 

Axial  momentum  transported  across  the  bottom 
boundary  of  a  cell. (See  Appendix  B) 

AMPY 

(PH2) 

g 

Mass  transported  across  top  boundary  of  a 
cell.  (See  Appendix  B) 

AMUR 

(PH2) 

cm-g/sec 

Radial  momentum  transported  across  right 
boundary  of  a  cell.  (See  Appendix  B)  > 

AMUT 

(PH2) 

-  - 

cm-g/sec 

Radial  monentum  transported  across  top 
boundary  of  a  cell.  (See  Appendix  B) 

.AMVR 

(PH2) 

-  - 

cm-g/scc 

Axial  momentum  transported  across  top 
boundary  of  a  cell.  (See  Appendix  B) 

AMX 

B.C. 

(KDX) 

& 

•Total  mass  in  a  cell.  (IMAX  by  JMAX  array) 

BBAR 

=  2(149) 

Used  in  CDT.  An  INPUT  parameter  used  in 
local  sound-speed  calculation  for  all 
materials  other  than  a  polytrcpic  gas. 

(Local  sound-speed  for  material  i  in  cell  K 
is  approximated  as  (C^; )  *  BBAR  *  *'P(k), 
where  the  coefficient  B  is  obtained  by 
determining  a  typical  slope  for  the  isen- 
tropes  in  Ref.  3  and  using  the  relation 

C  *  V /'-"dr  / «1V  to  evaluate  b  at  a  particular 
point. 

BBOUN’D 

=>2(74) 

"  • 

ergs 

•Calculated  in  PH5.  Printed  ir.  EDIT  under 
"Elastic  Plastic  Work."  Total  work  done  by 
the  clastic  and  plastic  stresses. 

BOTM 

-2(39) 

-- 

& 

Calculated  in  P112.  Printed  in  EDIT.  Total 

mass  transported  across  bottom  of  grid. 
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Variable 

Name 

Location 

Dimension 

Units 

Definition 

BOTMU 

BZ  (64) 

-- 

cm-g/sec 

Calculated  in  PII2,  Printed  in  EDIT.  Total 
radial 'momentum  transported  across  bottom  o: 
grid. 

BOTMV 

sZ  (40) 

*  * 

cm^g/sec 

Calculated  in  PH2.  Printed  in  EDIT.  ’  Total 
axial-momentum  transported  across  bottom  of 
grid. 

CNAUT 

MXCELL 

(30) 

cm/sec 

Used  in  CDT.  Approximate  sound-speed  of 
nineteen  materials  defined  in  a  DATA  state¬ 
ment  in  BLOCK. 

% 

/LScaPa. 

COi  Y  KllOZi  YAi/0Oi  * 

CVIS 

=Z(27) 

*  " 

"*  * 

INPUT  parameter.  Used  to  describe  the 
bottom  boundary -condition  in  PHI,  PH 2  and 
PH3.  Botton  boundary  is  transmittive  when 
CVIS  =  -1.,  reflective  when  CVIS  =  0. 

ere 

B.C. 

*  - 

*  * 

Defined  in  INFACE.  Used  in  FLUX  to  define 
flux  terns  of  cells  that  become  mixed  after 
first  subcyclc  of  INFACE, 

CYCLE 

“Z  (2) 

•  - 

-  - 

Used  in  INPUT,  SETUP,- CDT,  EDIT.  Cycle 
number  (an  integer  value  in  floating  point 
form) . 

CYCMX 

eZ  (69 ) 

INPUT  parameter.  Equals  number  of  passes 
through  INFACE  on  each  cycle.  ,  Suggested 
minimum  of  2  and  maximum  of  8.  Used  to 
minimize  transport  noise  near  interfaces. 

CYCPH3 

=Z(70) 

-- 

INPUT  parameter.  Number  of  times  to  sub¬ 
cycle  PH3 .  If  CYCPH5  <.0,  PH3  is  omitted. 

CZERO 

B.C. 

(KMX) 

dynes/ cm ^ 

Value  of  Yq  for  yield  strength  calculation. 
Defined  by  input  cards  for  each  material 
package  in  the  grid.  Used  in  STRNG. 

(See  STRONG) 

DDX 

B.C. 

(IDXX) 

cm 

An  array  equivalcnced  to  the  DX  array  such 
that  DDX(l)  =  DX(O) . 

DDY 

B.C. 

(JDXX) 

cm 

An  array  equivalenced  to  the  DY  array 
such  that  DDY(l)  =  DY(0). 

DELEB 

(PH2) 

_  - 

ergs 

Total  energy  associated  with  mass 
transported  across  botton  boundary  of 
a  cell.  (See  Appendix  B) 

DELET 

(PH2; 

■  * 

ergs 

Total  energy  associated  with  mass  trans¬ 
ported  across  top  boundary  of  cell.  (See 
Appendix  B) 

DELI 

(PIU,  P1I2) 

-- 

ergs/g 

Change  of  specific  internal  energy  of  a 
cell. 
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Variable 

Name 


Location  Dimension 


Units 


Definition 


DELM 

(PH2) 

•  - 

Z 

Total  change  in  mass  of  a  cell. 

DELU 

(PHI,  PH2) 

.. 

cm/sec 

Change  of  radial  velocity  of  a  cell. 

DELV 

(PHI,  PH2 ) 

-- 

cm/sec 

Change  of  axial  velocity  of  a  cell. 

DM  IN 

=2(24) 

INPUT  parameter.  Allowable  relative  error 
in  energy  sum.  When  relative  error  is 
>  DM IN ,  calculation  is  terminated.  The 
energy  error  is  computed  in  EDIT.  If  every¬ 
thing  is  wgrking  right  you  should  be  able 
to  use  lO'-*  for  DMIN. 

DPDR 

B.C. 

-- 

ergs/g 

Used  in  CDT  when  iterating  for  pressures. 
Computed  in  EQST  or  EOSXPL.  =(3P/3p)£ 

DT 

-Z(3) 

-• 

sec 

Time  step.  Calculated  in  CDT. 

DTFACT 

(PH3) 

•  - 

*  *■  ' 

Factor  used  in  calculating  a  variable  time 
step  when  subcycling  the  PHI  and  PH3 
calculatioa . 

DTFIX 

=Z (143) 

•  _ 

sec 

■  INPUT  parameter.  Used  in  CDT.  If 

DTFIX. GT. DTMIN,  then  DT=min(DT , DTFIX) . 

Used  to  run  with  a  constant  time  step. 

DTMIN 

“Z  (144) 

*  * 

sec 

INPUT  parameter.  Minimum  time  step;  used 
in  CDT.  After  STAB  =  FINAL,  if  DT  <  DTMIN 
execution  is  stopped. 

DTNA 

=2(263 

-- 

sec 

DT  from  previous  time  cycle.  Used  in 
.  INPUT,  CDT,  and  EDIT, 

DTSTR  ■ 

Local 

-- 

sec 

Computed  in  PHI.  Time  Step  for  PHI  and 

PH3  subcycles. 

DX 

B.C. 

(IDX) 

cm 

The  radial -dimension  of  cells.  Equi' 
valenced  to  DDX  such  that  DUX(l)  =  DX(0), 

DXF 

*Z (136) 

'  - 

cm 

An  INPUT  parameter  used  to  calculate  the 

DX  array  if  the  radial  dimension  of  the 
cells  is  uniform. 

DY 

B.C. 

(JDX) 

cm 

The  axial-dimension  of  cells.  Equi- 
valenced  to  DDY  so  that  DDY(l)  =  DY(0). 

DYF 

“Z (137) 

*  - 

cm 

An  INPUT  parameter  used  to  calculate  the 

DY  array  if  the  axial  dimension  of  the 
cells  is  uniform. 

EAMMP 

(PH2) 

--  . 

ergs/g 

Specific  internal  energy  of  mass  transported 
across  the  right  edge  of  cell. 

EAMPY 

(PH2) 

-- 

ergs/g 

Specific  internal  energy  of  mass  transported 
across  top  of  cell. 
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Variable 

Name 

Location 

Dimension 

Units 

Definition 

ECK 

-Z(76) 

Used  in  EDIT.  Relative  error  in  energy  sum: 
(E  Ek  -  ETII  )/ET!J; 

where  is  total  energy  in  cell  K.  If 

( ECK |  >  D.MIN,  execution  is  stopped. 

EMIN 

“Z (85} 

ergs/g 

INPUT  parameter.  Minimtm  specific  internal 
energy  of  an  ideal  gas  or  explosive  to  be 
used  in  the  pressure  iteration  in  CDT. 
Usually  EMIN  *  106. 

EMOS 

*Z(150) 

-- 

ergs 

Calculated  in  PH2.  Printed  in  EDIT.  Energy 
transported  across  bottom  of  mesh  . 

EMOR 

flZ (135) 

-  - 

ergs 

Calculated  in  PH2.  Printed  in  EDIT.  Energy 
transported  across  right  side  of  mesh. 

EMOT  ' 

=Z (146) 

-- 

ergs 

Calculated  in  PU2.  Printed  in  EDIT.  Energy 
transported  across  top  of  mesh. 

ENERGY 

B.C. 

ergs/g 

Defined  in  CDT  a;;  specific  internal  energy 
of  a  material.  Used  in  EQST  and  EOSXPL 
where  P  -  f (ENERGY, RHON) . 

EOB 

eZ (134) 

—  ** 

ergs 

Calculated  in  PHI.  Printed  in  EDIT.  Energy 
change  due  to  work  done  at  bottom  boundary 
of  grid. 

EOR 

“Z (132) 

ergs 

Calculated  in  PHI.  Printe'd  in  EDIT.  Energy 
change  due  to  work  done  at  right  boundary 
of  grid. 

EOT  . 

-ZC133) 

ergs 

Calculated  in  PHI.  Printed  in  EDIT.  Energy 
change  due  to  work  done  at  top  boundary  or 
grid. 

ERDUMP 

B.C. 

-* 

-- 

Used  in  EDIT  and  ERROR.  Flags  EDIT  to  do 
only  a  tape  dump  on  an  error  exit. 

ERR 

(PH3) 

(3,  JDX) 

1/sec 

Used  in  PH3 .  Radial  component  of  the 
instantaneous  strain-rate  deviators. 

ERZ 

(PH3) 

(3,  JDX) 

1/sec 

Used  in  PH3.  Shear  component  of  the 
instantaneous  strain-rate  deviators. 

ESA 

(EQST) 

(30) 

Defined  in  DATA  Statement.  Values  of  "a" 
in  equation  of  state  for  19  materials. 

[=(y-l)  when  using  perfect  gas  equation 
of  state. J 

ESALPH 

(EQST) 

(30) 

Defined  in  DATA  Statement.  Values  of  "a" 
in  equation  of  state  for  19  materials. 

ESB 

(EQST) 

(30) 

Defined  in  DATA  Statement.  Values  of  "b" 
in  equation  of  state  for  19  materials. 
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Variable 

Name 


Location  Dimension 


Units 


Definition 


ESBETA 

(EQST) 

(30) 

Defined  in  DATA  Statement.  Values  of  "8” 
in  equation  of  state  for  19  materials. 

ESCAPA 

(EQST) 

(30) 

-4ynes/cm2 

Defined  in  DATA  Statement.  Values  of  "A" 
in  equation  of  state  for  19  materials. 

ESCAPB 

CHQST) 

(30) 

dynes/cm2 

Defined  in  DATA  Statement.  Values  of  *'B" 
in  equation  of  state  for  19  materials. 

ESES 

(EQST) 

(30) 

ergs/g 

Defined  in  DATA  Statement.  Values  of  "E  " 
in  equation  of  state  for  19  materials. 

ESESP 

(eqst) 

(30) 

ergs/g 

Defined  in  DATA  Statement.  Values  of  "Eg" 
in  equation  of  state  for  19  materials. 

ESEZ 

(EQST) 

(30) 

ergs/g 

Defined  in  DATA  Statement.  Values  of  "Eg" 
in  equation  of  state  for  19  materials. 

ETH 

**2(13) 

ergs 

Theoretical  value  of  total  energy  in  the 
mesh.  Calculated  in  SETUP  initially; 
updated  in  APLYPR,  PHI,  PH2  and  R.MVMAT. 

EVAPEN 

°z  (10 1) 

ergs 

Calculated  in  PH2  and  CDT .  Printed  in 

EDIT.  Sum  of  energy  lost  through  "evapora¬ 
tion"  of  mass  left  in  cells  due  to  round-off 
error.  Adjusted  in  CDT  when  "evaporating" 
isolated  cells.  Initialized  in  SETUP. 

EYAPM 

«Z(100) 

*  - 

E 

Calculated  in  NI2  and  CDT.  Printed  in  EDIT. 
Sum  of  mass  lost  through  "evaporation". 

See  EVAPEN. 

EVAPMU 

-Z(102) 

*  * 

cm-g/sec 

Calculated  in  PI 1 2  and  CDT.  Printed  in  EDIT. 
Sum  of  radial  momenta  lost  through 
.  "evaporation" .  See  EVAPEN. 

EVAPMV 

-Z(103) 

cm-g/sec 

Calculated  in  PH2  and  CDT.  Sum  of  axial 
momenta  lost  through  "evaporation".  See 
EVAPEN. 

EZPH2 

■Z  (104) 

■  - 

ergs/g 

Surfi  of  specific  internal  energy  fluxes  that 
are  less  than  IJMIN-  and  are  set  to  zero  in 

PH2.  Printed  in  EDIT. 

EZZ 

(PH3) 

(3,  JDX) 

1/sec 

Used  in  PH3.  Axial  component  of  the  instan¬ 
taneous  strain-rate  deviators. 

FINAL  • 

*»Z  (113) 

•  * 

-  - 

INPUT  parameter.  Maximum  value  of  stability 
fraction  (STA3) .  If  FINAL  -  0;  the  stability 
fraction  will  be  constant.  Used  in  CDT. 

FLEFT 

(PH2) 

(JDX) 

cm-g/sec 

Radial  momentum  of  mass  transported  across 
left  side  of  cell.  Equivalenced  to  UL  array. 
(See  Appendix  B) 

FRACRT 

B.C. 

(MDKMX) 

cm2 

The  fractional  area  of  the  right  cell  boundary 
associated  with  a  given  material  in  a  mixed 
cell.  Used  to  compute  mass  flux  of  that 
material. 
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Variable 

Name 


Location  Dimension 


Units 


Definition 


FRACTP 

B.C. 

(MD.VMX) 

cm2 

The  fractional  area  of  the  top  cell  boundary 
associated  with  a  given  material  in  a  mixed 
cell.  Used  to  compute  mass  flux  of  that 
. material . 

CAMC 

(PH2) 

(JDX) 

g 

Mass  transported  across  left  side  of  cell. 
Equivalenccd  to  PL  array.  (See  Appendix  8) 

GAMMA 

-2(62) 

Defined  during  SETUP.  GAMMA  *  ESA+1  in 
ideal  gas  or  explosive  equation  of  state. 
Used  in  PHI  to  partition  energy  in  mixed 
cells. 

I 

°Z(88) 

■ 

-- 

Used  in  most  subroutines  as  index  in  radial 
direction . 

ICSTOP 

-2(7) 

*  - 

-  - 

INPUT  parameter.  Used  in  EDIT.  Execution 
stops  on  ICSTOP  cycle  when  stopping  on  a 
specified  cycle  rather  than  time. 

IGM 

*2(21) 

-  - 

-  - 

INPUT  parameter.  If  IGM*Q ,  problems  are 
run  using  cylindrical  geometry,  if  IGM-1, 
rectangular  geometry. 

1MAX 

®Z  (33) 

-  -  • 

-- 

INPUT  parameter.  .Number  of  columns  in 
mesh . 

IMAXA 

*Z(34)  . 

IHAX  ♦  1. 

INTER 

3Z (37) 

INPUT  parameter.  On  print  cycle  if 

INTER  f  0,  EDIT  will  print  after  CDT 
and  PHI.  If  INTER  =  7,  energy  totals 
are  printed  in  PH2  in  addition  to  the  extra 
EpiT  prints.  . 

IPCY.CL  ... 

-2(49) 

INPUT  parameter.  Used  in  EDIT.  The  number 
of  cycles  between  EDIT  prints  when  editing 
on  cycles  rather  than  time. 

IPR 

»z  as) 

•  -  - 

-  - 

Maximum  number  of  iterations  to  be  performed 
by  CDT  to  achieve  pressure  equilibration 
between  materials  in  nixed  cells. 

IVARDX 

“Z (83) 

-- 

--  • 

When  IVARDX  »  1,  radial  dimension  of  cells 
is  variable. 

1VARDY 

-2(54) 

-- 

-- 

When  IVARDY  *  1,  the  axial  dimension  of 
cells  is  variable. 

n 

•2(47) 

-  - 

-- 

11  is  initialised  in  SETUP  and  is  used  to 

limit  calculation  in  the  radial  direction  to 
an  "active  mesh",  Beyond  II  nothing  is  vet 
disturbed  from  tlu:  initial  conditions.  II 
is  specified  initially  as  (2  ♦  the  column- 
number  of  the  last  column  in  which  there  is 
a  non-zero  velocity  or  internal  energy). 

11  is  increased  automatically  as  inactive 
cells  become  active.  However,  II  is  never 
larger  than  I MAX. 


i 


J 


n 

u 

.€> 


tr- 
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Variable 

Name 

Location 

Dimension 

Units 

Definition 

12 

-2(48) 

-• 

-- 

Initialized  in  SETUP.  Like  II  but  for  axial 
disturbance  limit. 

13 

B.C . 

*  * 

Used  in  EDIT  as  a  flag  for  "short"  or  "long" 
prints.  13  !  number  of  columns  of  grid  whos< 
U,  V,  P,  AIX,  etc.  are  listed. 

J 

°Z  (89) 

-- 

-- 

Used  as  row-index  in  most  subroutines. 

JMAX 

-Z(3S) 

-- 

-- 

INPUT  parameter.  Number  of  rows  in  the 
grid. 

JMAXA 

=  Z(36) 

-- 

-- 

Calculated  in  SETUP,  (JMAX  ♦  1). 

K 

-2(90) 

-- 

-- 

Used  as  cell-index  in  all  subroutines. 

K  =  (J-1)*1MAX  ♦  I  ♦  1. 

KA 

B.C. 

--  • 

-- 

Used  as  cell-index  for  cell,  above  cell  K. 

KB 

B.C. 

-- 

-- 

Used  as  cell-index  for  cell  below  cell  K. 

KL 

B.C. 

-- 

-- 

Used  as  cell-index  for  cell  to  left  of 
cell  K. 

KMAX 

»Z(37) 

-- 

-- 

Calculated  in  SETUP  (= I MAX* JMAX  ♦  1). 

K-index  of  last  cell  in  grid. 

KMAXA 

-2(38) 

-- 

Calculated  in  SETUP  (KMAX  +  1). 

KR 

B.C. 

.  -- 

-- 

Used  as  cell  index  for  cell  to  right  of 
cell  K. 

KSPACE 

(EDIT) 

-  -  . 

— 

Used  for  spacing  printed  output. 

KUNITR 

=  2  (14) 

-- 

-- 

Unit  number  of  tape  read  by  INPUT. 

KUNITW 

“2(17) 

-- 

-- 

Unit  number  of  tape  written  on  by  SETUP  and 
EDIT. 

M 

“Z (91) 

r- 

Usually  used  as  an  index  for  cell  K  in  the 
mixed  cell  arrays.  M--  MFLAG(K) -100 . 

NA 

B.C. 

-- 

-- 

Usually  value  of  MFLAG  for  cell  above 
cell  K. 

MAPS 

*2(4  2) 

INPUT  parameter  that  determines  the  printing 
of  symbolic  maps  on  EDIT  cycles.  If  MAPS  -  i 
maps  are  not  printed;  if  MAPS  -  1,  they  are 
printed. 

MAT 

B.C. 

(30) 

Material  code  number  for  the  liner. 

MAT(l)  a  3  indicates  th“c  material  in  package 
one  is  iron.  See  list  in  Section  2.5  of 
this  report  or  comments  in  listing  of  EQST 
subroutine. 

MATXPL 

n 

r-o 

i— 

o 

O' 

*  * 

•  • 

If  MAT(2)*20,  package  two  is  an  ideal  gas 
(explosive).  MATXPL  thou  identifies  the 
explosive . 

MFK 

B.C. 

-- 

-- 

Usually  -  Ml:LAC(K). 
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Variable 

Name 


Location  Dimension 


Units 


Definition 


MFLAC 

B.C. 

(KDX) 

Flag  for  each  cell  that  indicates  whether  it 
is  pure  or  mixed.  If  MFLAG(K)  <  100  cell  is 
pure.  If'MFLAG(K)  >  100  cell  is  mixed  and 
flag  gives  storage  location  in  mixed  cell 
arrays  for  cell  k  (see  Appendix  C): 

MO 

B.C. 

Defined  in  CALFRC.  MFLAG(K)  of  cell  K  before 
it  becomes  mixed.  Used  in  NEIVMIX  to  define 
mixed  cell  variables  for  cell  K. 

MR 

B.C. 

Value  of  MFLAG  for  cell  on  right. 

N 

B.C. 

-- 

-- 

In  EQST,  N  is  material  code  number  trans¬ 
ferred  from  CDT. 

N6 

‘  =  Z  (56) 

*  “ 

INPUT  parameter.  Used  in  CDT.  Negative 
pressures  are  allowed  in  cells  above 

J  =  N6. 

N10 

=2(60) 

--  ‘ 

-- 

Used  in  CDT  to  identify  I-index  of  cell 
which  controls  the  time  step. 

Nil 

=  Z(61) 

-- 

Used  in  CDT  to  identify  J-index  of  cell 
which  controls  the  time  step. 

NAPRND 

“Z  (23)  ‘ 

Used  in  APLYPR  to  identify  the  last  free 
surface  tracer  describing  the  section  of 
the  free  surface  which  will  have  3  pressure 
force  applied  to  it. 

NAPRST 

=2(23) 

.  “  ” 

*  * 

Similar  to  NAPRST  except  it  identifies  the 
first  tracer.  NAPRST  =  0  means  that  no 
pressure  force  is  to  be  applied. 

NC 

=Z(30) 

-- 

Cycle  number.  Set  initially  to  -1  in  input. 
Incremented  thereafter  in  CDT. 

-N  DUMP  7 

=  2(6) 

INPUT  parameter.  Used  in  EDIT  to  control 
frequency  of  tape  dumps.  e.g.,  tape  dump 
will  occur  on  every  EDIT  print  when  NDUMP7 
=1,  or  on  every  fifth  EDIT  print  when 

NDUMP7  =5. 

NECYCL 

=  2(77) 

*  - 

Defined  and  printed  in  EDIT.  The  cycle 
during  which  the  largest  relative  error  in  the 
energy  sum  was  computed. 

NERR 

B.C. 

m  * 

*  * 

Used  in  ERROR  as  exit  flag.  Prevents  ERROR 
from  being  called  more  than  once  during  a 
single  run. 

NFRELP 

.  “2  (S) 

INPUT  parameter.  Used  in  EDIT  to  control 
frequency  of  "long"  prints,  A  "long"  print 
will  occur  with  every  UIUT  if  NFRELP  =  1; 
with  every  fifth  EDIT  i£  NFRELP  *  5  (see  13). 

NK 

B.C.  . 

-- 

-- 

Tells  which  statement  number  of  a  subroutine 
is  near  the  C3ll  to  ERROR. 

Variable 

Name  Location 

Dimension. 

Units 

Definition 

NKDX4  °Z(117) 

-- 

-- 

Number 
are  in 

of  interior  tracers  (XP,  YP)  that 
the  problem. 

NMAT 

eZ  (6  8) 

•  • 

•  • 

INPUT  parameter.  Equals  number  of  material 
packages  in  the  problem. 

NMP 

B.C  . 

(N’MXX) 

-- 

NMP(N)  is  the  number  of  interface  tracers 
describing  material  package  N. 

NMXCLS- 

=Z(73) 

• 

INPUT  parameter.  Equals  the  maximum  number 
of  mixed  cells  that  SETUP  or  N'EKMIX  can 
generate.  This  number  should  coincide 
with  the  dimension  of  variables  in  the 

MXCELL  common  block. 

NODUMP  ' 

=  2(96) 

• 

INPUT  parameter.  Used  in  EDIT.  When 

NODUMP  =  1 ,  no  tape  dumps  arc  made  except 
by  SETUP  on  cycle  0.  Allows  user  to  re¬ 
start  a  problem  without  writing  on  the 
restart  tape. 

NPRINT 

B.C. 

-- 

-- 

NPRINT  =  1  during  the  cycle  on  which  EDIT 
prints  and  checks  the  energy  error. 

NR 

.  B.C. 

-- 

-- 

Identifies  which  subroutine  calls  ERROR. 

Used  in  ERROR  for  printing  error  message. 

NTCC 

-Z(S1) 

--  . 

-- 

INPUT  parameter.  When  NTCC  =  1  passive 
cell  centered  tracers  arc  used. 

NTCRSV 

“Z  (97) 

"  • 

Defined  in  SETUP.-  Used  in  MOVTCR  to  iden¬ 
tify  the  tracer  at  the  lower  left  hand  corner 
of  the  liner. 

NTPMX 

“2 (78) 

INPUT  parameter.  Equals  the  maximum  number 
of  tracers  that  SETUP  will  generate  to  cir¬ 
cumscribe  a  single  material  package.  The 
value  of  this  number  should  not  exceed  the 
dimension  of  the  TX  and  TY  arrays. 

NUMSCA 

-2(43) 

*  " 

•  “ 

INPUT  parameter-  Number  of  times  the  print 
interval  is  to  be  rescaled.  Used  in  EDIT. 

See  PRUELT  for  further  details. 

NUMSP 

“2(4) 

-* 

Used  in  EDIT  to  count  the  number  of  prints 
(short  or  long)  since  the  last  tape  dump. 

NVOID 

B.C. 

*  * 

®  • 

Defined  in  INPUT:  =NMAT+1,  number  of 
material  packages  +  1.  The  number  of  the 
void  package  if  there  is  one. 

P 

B.C. 

(KDXP) 

dynes /cm^ 

Cell  pressure.  (iMAX  by  JMAX  array). 

Calculated  in  EQST.  UsoJ  by  Pill.  (NOTE: 
the  P-storagc  space  is  used  for  other 
information  in  INTACT  and  ITU . J 

PIDY 

-Z(8) 

-- 

-- 

-  n (3.1415927) 

PL 

B.C. 

(JDXP) 

dynes/cm^ 

Used  in  PHI  for  saving  pressures  on  left 
side  of  cell.  (.Storage  used  for  other 

purposes  in  CUT  and  l‘U2.)  s 
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Variable 

Name 


Location  Dimension  Units 


Definition 


PLW 

MXCELL 

(NMX) 

ergs 

An  array  that  stores  the  total  work  done  by 
plastic  stresses  in  each  material  package. 
Calculated  in  1*113 ;  printed  in  EDIT. 

PM  IN' 

aZ  (86) 

-- 

dynes/cm2 

Used  as  a  pressure  cut-off.  If  the  |P| 
is  less  than  PM  IN',  P  =  0. 

PRCNT 

-2(16) 

Convergence  requirement  for  equilibrating 
pressures  in  a  mixed  cell.  If 
(?.  -  F  v 

—  j  PRCNT  for  all  materials  (i)  in 
in  K,  materials  in  cell  arc  considered  to 

be  pressure  equilibrated  and  P(K)  =  P. 

PRDELT  =Z(4S)  sec  INPUT  parameter.  Gives  the  initial  time 

interval  between  EDIT  prints.  There  are 
five  parameters  which  control  printing 
frequency:  PRDELT,  IPCYCL,  PRLIM ,  PRFACT, 
and  NUMSCA.  If  the  user  is  printing  on  time 
(PRDELT  f  0  and  ITCYCL  =  0),  DT  will  be 
adjusted  so  that  a  print  will  occur  exactly 
every  PRDELT  seconds.  If  the  user  is 
printing  on  cycles  (PRDELT  =  0,  IPCYCL  /  0) » 
a  print  will  occur  every  IPCYCL  cycles. 
PRLIM,  PRFACT  and  NUMSCA  arc  used  to  in¬ 
crease  the  print  interval.  PRLIM  is  the 
time  (or  cycle)  at  which  PRDELT  (or  PICYCL) 
and  PRLIM  are  multiplied  by  PR FACT.  The 
new  value  of  PRLIM  establishes  the  next 
time  (or  cycle)  when  the  print  interval 
will  be  rescaled.  This  process  continues 
at  most  NUMSCA  times.. 

EXAMPLE:  You  wish  to  print_cvcry  1  *  10"® 
sec  until  you  reach  I  *  10"'  sec,  then  every 
1  x  10"'  see  until  1  x  10"6  sec  and  every 
1  x  10 sec  thereafter: 

PRDELT  «  1.  x  10'? 

PRLIM  =  i.  x  10 ”7 

PRFACT  -  10. 

NUMSCA  *  2. 


PRESUR 

B.C. 

dynes/ cm2 

Defined  in  EQST :  pressure  =  f(o,  E) .  Used 
in  CDT  to  define  P(K)  in  the  case  of  pure 
cells,  and  in  the  case  of  nixed  cells  to 
define  the  pressure  of  a  material. 

PRFACT 

”2(46) 

INPUT  parameter.  Used  in  EDIT  as  a  factor 
for  rescaling  PRUELT ,  IPCYCL  and  PRLIM  when 
PRLIM- fine  or  cycle  is  reached  (see  PRUELT) 
Should  be  >  1 . 

PRLIM 

*Z (4  4) 

■  - 

-  - 

INPUT  parameter.  Tine  or  cycle  at  which 
to  rescale  PKDHLT  (or  IPCYCL)  and  PRLIM  by 
PRFACT  (see  PR PEL I) .  / 

PROB 

=  2(1) 

-- 

-- 

INPUT  parameter.  Identifying  pro!)  1  era 
number.  Should  be  between  0  and  IU0. 

PRTLME 

*2(131) 

-• 

sec 

Initially  set  to  I'KDJ.LT  in  INPUT.  There¬ 
after  calculated  in  EDIT.  Mien  T  ■  PUT I ME, 
it  is  time  to  print. 
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Variable 

Name 

Location 

Dimension 

Units 

Definition 

RATIO 

(CDT) 

-- 

sec 

Used  in  calculation  of  DT:  «  DX  •  DY/ 

( DX • I V |  ♦  DY-|U|). 

RELERR 

(EDIT) 

Used  for  storing  and  printing  maximum 
relative  error  in  the  energy  sum. 

R1I0 

MX CELL 

(NMXX  ,MD.X) 

S/cm3 

The  density  of  materials  in  mixed  cells. 
Updated  in  CDT  and  PH2.  If  cell  K  is  mixed 
and  M  =  MFLAG (K)  - 100  ,  then  RIIO(l,M)  is  the 
density  of  material  1  in  cell  K.  (See 
Appendix  C) 

RHOIN 

B.  C« 

(NMX) 

g/cm3 

The  initial  densities  of  material  in  the 
packages.  Defined  by  input  cards  read  by 
SETUP. 

RHOW 

B  .C . 

-  • 

g/cm3 

Density  of  material.  Defined  in  CDT  and 
used  in  EQST  to  define  pressure:  P  « 
f  (ENERGY  , RM05V)  . 

RHOZ 

MXCELL 

(30) 

g/cm3 

Defined  in  DATA  statement  in  BLOCK.  Normal 
density  for  19  materials. 

RMU 

B.C. 

(NMX) 

dynes/cm2 

Rigidity  modulus  values  for  material 
packages.  Defined  by  input  cards  read  by 
SETUP.  Used  in  PI13. 

ROEPS 

“2 (110) 

«  _  — 

INPUT  parameter.  A  round-off  epsilon  used 
in  calculating  cutoffs  of  energy*  velocity* 
pressure  and  mass  flus. 

RTIPI 

=2(28) 

-- 

cm 

Radius  of  tip  of  liner.  (See  YCEXT) 

If  RTIPI  <  0,  liner  has  a  pointed  tip. 

RTM 

°2(S7) 

S 

Calculated  in  PH2.  Printed  in  EDIT.  Total 
mass  lost  out  right  side  of  grid. 

RTMU 

-2(10) 

-- 

ca-g/sec 

Calculated  in  PH2 .  Printed  in  EDIT.  Total 
radial -momentum  lost  out  right  side  of  grid. 

RTMV 

*2(58) 

-- 

cra-g/sec 

Calculated  in  PH2 .  Printed  in  EDIT.  Total 
axial -momentum  lost  out  right  side  of  grid. 

SAMMP 

MXCELL 

(NMX ,MDX) 

g 

The  flux  of  materials  across  right  boundary 
of  mixed  cells.  Calculated  in  l-'LUX;  used 
in  PH2  and  INFACE. 

SAMPY 

MXCELL 

(NMX.MDX) 

g 

The  flux  of  materials  across  the  top  boundary 
of  mixed  cells.  Calculated  in  FLUX;  used 
in  Pll2  and  INFACE. 

SDT 

B.C. 

-- 

sec 

Defined  in  INFACE.:  -  UT/CYCMX.  Tine  step 
for  each  subcycle  of  INFACE. 

SIE 

MXCELL 

(NMX.MDX) 

ergs/g 

The  specific  internal  energy  of  materials. 

in  nixed  cells.  !f  cell  K  is  nixed  and 
M  =  MFl.AG(K)- 100  ,  then  is  the 

specific  internal  energy  oi'  material  one 
in  cell  K  (see  XMASS,  UliO) . 
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Variable 

Name 

Location 

Dinens ion 

Units 

Definition 

SIENEK 

CPH2) 

•- 

crgs/g 

New  value  of  specific  internal  energy. 

SIGC 

B.C. 

um 

ergs 

Used  in  PH2  for  internal  energy  of  mass 
transported  across  left  side  of  cell. 

(Sec  Appendix  B) 

SIGMU 

mi) 

-- 

cm-g/sec 

Total  change  in  radial  momentum  of  a  cell. 

SIGMV 

(PH2) 

-- 

cm-g/sec 

Total  change  in  axial  momentum  of  a  cell. 

SLOP LX 

=2(31} 

-- 

-* 

Slope  of  inside  surface  of  liner  (see  YTOPI). 

SOLID 

(STRXG) 

~  • 

g/cm3 

=  RH0Zl  *  AMDMj..  In  PH3 ,  if  RHO*  <  SOLID 
for  one  material,  i,  in  cell  K,  stresses 
of  cell  K  are  set  to  zero. 

SRATIO 

(CDT) 

*- 

sec 

Used  to  calculate  DT.  Smallest  value  of 

RATIO  found  in  the  active  grid. 

SRR 

(PH3) 

(3,  JDX) 

dynes/cm2 

Temporary  storage  for  ceil  centered 
deviatoric  normal  stress  in  the  radial 
direction.  See  STRSRR. 

SRZ 

(PH3) 

(3  »■  JDX) 

dynes/cn2 

Similar  to  SRR  but  for  shear  stress.  See 
STRSK2 . 

SSI  EX 

B.C. 

(NMX) 

ergs/g 

Initial  specific  internal  energy  of  each 
package.  Defined  by  input  cards  read  in 
SETUP. 

STAB 

=2(139) 

* 

INPUT  parameter.  Used  in  CDT.  Initial 
value  of'stahi lit;/  fraction"  for  the  cal¬ 
culation  of  DT.  If  FINAL  =  0.,  STAB  is 
constant.  Otherwise  its  value  progresses 
from  STAB  to  FINAL  in  a  geometric  progression 
{'NOTE:  DT  =  STAB  *  SRATIO] 

STEZ 

B.C. 

(KMX) 

ergs/g 

Eq  for  each  material  package.  Used  in 
yield-strength  calculation  in  STRXG. 

Defined  by  Input  cards  read  by  SETUP. 

(See  STRENG) 

STk'l 

B.C. 

(NMX) 

dynes/cm2 

Y.  for  each  material  package.  Used  in 
yield-strength  calculation  in  STRXG. 

Defined  by  input  cards  read  by  SETUP. 

(See  STRENG) 

STRENG 

(STRXG) 

-- 

2 

dynes/cm 

Yield  strength  of  material: 

STRENG  *  (Y0  ♦  YjPJCl  -  E/Eg) 

y  3  p/pn  -  1.  If  STRENG  <_  0 . ,  stresses 
are  set-  to  0.  If  E  >  Kg,  STRENG  =  0  and 
stresses  arc  set  to  0. 


1 

i 

3 
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Variable 

Name 


Location  Dimension 


Units 


Definition 


STRENG  (continued)  .  .  YQ  is  CZERO, 

Yj  is  STK1 , 

P  is  pressure  of  material, 

E  is  specific  internal  energy  of  material, 
Eq  is  STEZ, 

p  is  density  of  material. 


Pq  is  normal  density  of  material. 


STRSRR 

ELPL 

(KDX) 

2 

dynes/cm 

The  cell  centered  deviatoric  normal  stress 
in  the  radial  direction  (I MAX  by  JMAX  array). 

STRSRZ 

ELPL  ' 

(KDX) 

2 

.  dynes/cm 

The  cell  centered  deviatoric  shear  stress  . 

(I MAX  by  JMAX  array). 

STRSZZ 

ELPL 

(KDX) 

dynes/cm2 

The  cell  centered  deviatoric  normal  stress 
in  the  axial  direction  (IMAX  by  JMAX  array). 

SUM 

B.C. 

-- 

-- 

Used  in  most  routines  as  working  storage 
when  summing  quantities. 

SUME 

CPU2) 

ergs 

Used  in  PH2.  Sums  energy  fluxes  ignored 
during  one  cycle.  Used  to  adjust  ETH, 
the  theoretical  energy  total. 

SZZ 

CPH3) 

(3 ,  JDX) 

2 

dynes/cm 

Similar  to  SRR  but  for  axial  direction. 

See  STRSZZ. 

T 

=  Z(84) 

sec 

Time.  Initially  defined'  in  INPUT.  Incre¬ 
mented  in  CUT.  Adjusted  in  EDIT  for 
printing  and  stopping  at  specified  time 
values. 

TAU 

B.C. 

(IDX) 

cn2 

Initially  defined  in  SETUP.  Area  of 
'cell  face: 

“  m[X(I)2  -  X(I-1)21  for  cylindrical  problems 

"  X(I)  -  X(I-l)  for  rectangular  problems. 

Used  in  most  subroutines.  Sec  IGM. 

TAUDTS 

(PHI) 

-- 

cm2/ sec 

=  TAU  *  DT. 

TH03 

(PH3) 

-- 

1/sec 

°  1/5(ur  *  vz  *  ?)• 

TJMOB 

=  2(98) 

g 

Total  jet  mass  which  has  passed  out  the 
’  bottom  of  the  grid.  Computed  and  printed 
in  PJI2  every  cycle. 

TOPM 

«Z(63) 

-  - 

g 

Calculated  in  PH2.  Printed  in  EDIT.  Total 
mass  transported  across  top  of  grid. 

TOPMU 

=  2(9) 

-- 

cm-g/sec 

Calculated  in  PII2.  Printed  in  EDIT.  Total 

radial -momentum  transported  across  top  of 
grid. 
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Variable 

Kane 

Location 

Dimension 

Units 

Definition 

TOPMV 

CZ  (66) 

— 

cm-g/scc 

Calculated  in  PII2.  Printed  in  EDIT.  Total 
3xial -momentum  transported  across  top  of 
grid. 

TRIAL 

(CDT) 

-- 

cm/sec 

Used  in  CDT.  Maximum  sura  in  the  entire  grid 
of  a  cell's  sound-speed  and  largest  velocity- 
component.  Printed  in  CDT  as  MAX CUV . 

TSTOP 

aZ(5Q) 

•  * 

sec 

INPUT  parameter.  Value  of  T  at  which 
execution  stops  when  stopping  on  time  rather 
than  cycles. 

TKOPI 

B.C, 

— 

Calculated  in  INPUT.  =  2n. 

TX 

B.C. 

(NMXX.NTPX) 

cell  . 

X-cocrdinatcs  of  tracer  particles  circum¬ 
scribing  the  material,  packages. 

TY 

B.C. 

(NMXX,XTPX) 

cell 

Y-coardinates  of  tracer  particles  circum¬ 
scribing  the  material  packages. 

U 

B.C. 

(KDX) 

cm/scc 

Radial  velocity  of  cell  (I MAX  by  JMAX  array), 

UAMMP 

(PH2) 

-- 

cm/ sec 

Radial  velocity  of  mass  transported  across 
right  cell  edge. 

UAMPY 

(PH2) 

-- 

cm/sec 

Radial  velocity  of  mass  transported  across 
top  cell  edge. 

UEFF 

(MOVTCR) 

cm/sec 

Radial  velocity  of  a  tracer  particle,  com¬ 
puted  from  an  area-weighted  averaging  scheme. 

UL 

(PHI) 

(JDX2) 

cm/sec 

Radial  velocity  at  left  boundary  of  cell 

K.  UL  array  equivalcnccd  to  the  FLEFT  and 
YMAC  arrays  used  in  PH2. 

UMIN 

=2(129) 

-- 

cn/sec 

Calculated  in  CDT.  Used  as  velocity  cutoff 
in  PH2 :  =  (RJcPo)  4  (souno- speed  "  velocity) E 

UK  xxx 

=2(xxx) 

-- 

-- 

Unused  2 -storage 

UXEW 

(PH2) 

-- 

cn/sec 

New  radial  velocity  of  cell  K. 

uox 

(PH3) 

-- 

cm/sec 

*  2U  /(Xj  ♦  Xj.j).  Used  to  define  strain 
rates  of  cell  K. 

URR 

B.C. 

-- 

cm/sec 

Used  in  FLUX  and  PH2.  Terrorary  storage 
for  radial  transport  velocity  at  right 
boundary  of  cell  K. 

UUR 

B.C. 

(NMX) 

cn/sec 

Initial  radial  velocities  of  packages. 
Defined  by  input  cards  read  in  SL7U?. 

UVMAX 

=Z(S0) 

-- 

cm/sec 

Calculated  in  CDT.  Maximum  velocity 
in  grid. 

V 

B.C. 

(.KDX) 

cn/sec 

Axial  velocity  of  cell  (I MAX  by  JMIX  array). 

Variable 

Sane 

Location  Pir.ension 

Units 

Definition 

V ABOVE 

B.C. 

-- 

cin/sec 

Used  in  FLUX  and  PH2.  Temporary  storage  for 
axial  transport  velocity  at  top  boundary  of 
cell. 

VAMMP 

(PH2) 

-- 

co/sec 

Axial  velocity  of  mass  transported  across 
right  edge  of  cell. 

VAMPY 

(PJI2) 

*- 

cm/sec 

Axial  velocity  of  nass  moving  across  top 
edge  of  cell. 

VBLO 

(PHI) 

-- 

cm/sec 

Axial  velocity  at  bottom  boundary  of  cell. 

VEFF 

(MOVTCR) 

cra/sec 

Axial  velocity  of  a  tracer  particle,  computed 
from  an  area-weighted  averaging  scheme. 

VNEW 

CPH2) 

-- 

cffl/sec 

New  axial  velocity  of  cell  K. 

VOLCOF 

(PHI) 

(KMX) 

g2/sec2-cin4 

Used  in  PHI*  to'  compute  the  distribution  of 
internal  energy  tc  material  in  mixed  cells. 

VOLMAT 

MXCEIL 

(NMX ,MDX) 

c»3 

Volume  of  material  in  nixed  cells. 

»  XMASS (N ,M) /RHO(N #M) 

VOLP.T 

MXCEIL 

(NMX.MDX) 

o3 

Volume  flux  of  material  at  right  boundary 
of  cell.  (See  Appendix  B)- 

VOLTP 

MXCELL 

(KMX ,MDX) 

cm5 

Volume  flux  at  top  boundary  of  cell 
(See  Appendix  B) 

VOW 

(EQST) 

-- 

— 

•Pq/p. 

WA 

B.C. 

(KMX) 

cn/sec 

Initial  axial  velocities  of  packages. 

Defined  by  input  cards  read  by  SETUP. 

WFLAGF 

-Z(S1) 

Used  in  INPUT  and  EDIT.  Set  =  1.  on  first 
cycle  of  run  in  INPUT,  Triggers  an  EDIT 
print  on  first  cycle  of  every  run.  Reset 
to  0.  at  end  of  EDIT. 

WFLAGL 

=2(52) 

*  * 

— ’  * 

Used  in  MAIN  and  EDIT.  Flags  last  cycle. 

Set  *  1 .  in  EDIT.  Signals  BRLSC  to  call 
exit. 

KPB 

WPL 

WPR 

KPT 

(PHI) 

B.C. 

(PHI) 

(Pill) 

(JDX)  | 

dynes/cm2 

Used  in  PHI  to  store  the  interface  pressures 
at  cell  edges,  i.e.,  V.'PR  is  the  interface 
pressure  between  cell  K  and  cell  KR  (3  K*l) 

h’STB 

WSTL 

KSTR 

WSTT 

(PHI) 

B.C. 

(PHI) 

(PHI) 

(JDX)  | 

dynes/cm2 

Used  in  PHI  to  store  the  deviator  shear 
stress  at  each  interface  of  a  cell. 

ws 

h'SA 

KSB 

WSC 

WSX 

KSY 

B.C.  'I 
B.C.  j 
B.C.  1 
B.C.  | 
B.C. 
B.C.  ) 

1 

> 

-- 

Used  in  most  subroutines  for  local  working 
storage . 
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Variable 


Name 

Location 

Dir.ens  ion 

'  Units 

Definition 

KUB 

WUL 

WUR 

KUT 

(PHI.PH3) 

'  B  C 

(Pi’ll  *PII3) 
(PII1.PH3) 

(Ji»x)  | 

cm/sec 

Used  in  PHI  and  PII3  to  store  the  radial 
components  of  velocity  at  each  cell 
interface . 

WVB 

KVL 

KVR 

WVT 

(PHI ,PH3) 

B  C 

(PHL*PH3) 

(PH1.PH3) 

(JDX)  | 

cn/sec 

Used  in  Pill  and  PH3  to  store  the  axial 
components  of  velocity  at  each  cell 
interface. 

X 

B.C. 

(IDX) 

cm 

Distance  from  axis  to  outside  of  cell. 
Equivalcnced  to  XX  array  such  that  X(0) 

*  XX(1). 

XCONI 

=Z(99) 

era 

The  X-coordinatc  of  the  point  at  which  the 
curved  inside  tip  of  the  liner  meets  the 
straight  inside  section  of  the  liner. 

XIENRG 

-Z(140) 

*  * 

ergs 

Total  internal  energy.  Calculated  in  EDIT 
and  used  for  printing  labels  on  tracer 
particle  plots. 

XKENRG 

—  Z  (141) 

" 

ergs 

Total  kinetic  energy.  Calculated  in  EDIT 
and  used  for  printing  labels  on  tracer 
particle  plots. 

XMASS 

MXCELL 

(NMX.MDX) 

g 

The  mass  of  materials  in  mixed  cells 
(sec  SIE ,  RHO) . 

XMAX 

=Z(18) 

-- 

•  cm 

-  X (IMAX) . 

XP 

TRACRS 

(KDX4) 

cell 

X-coordinates  of  passive  cell -centered 
tracer  particles . 

XPLCNT 

=  Z ( 10  8) 

-- 

era 

The  axial  coordinate  of  the  point  where 
the  detonation  was  initiated. 

XPLDET 

(EOSXPL)  ' 

(18) 

era/ sec 

Defined  in  DATA  Statement.  Gives  detonation 
Velocities  for  IS  explosives. 

XPLEN'G 

(EOSXPL) 

(18) 

ergs/g 

Defined  in  DATA  Statement.  Values  of 
specific  detonation  energy  for  IS  explosives 

XPLGAM 

(EOSXPL) 

(18) 

•Defined  in  DATA  Statement.  Values  of  ”yM 
for  18  explosives. 

XPLRAD 

=  Z  (107) 

-  - 

cm 

The  initial  (T=0)  radius  of  the  detonation 
front.  If  a  plane  detonation  front  is  used, 

*  0. 

XPLRIIO 

(EOSXPL) 

(13) 

g/cm3 

Defined  in  DATA  Statement.  Values  of 
initial  densities  for  13  explosives. 

XPLVEL 

*2(109) 

-- 

cra/sec 

Velocity  at  which  explosive  package  is  being 

moved  through  grid.  Set  *  V \ A ( I ) . 
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Vari able 
Mane 

Location 

Dimension 

Uni  ts 

Definition 

XTENRG 

“2(142) 

-- 

ergs 

Total  energy.  Calculated  in  EDIT  and  used 
for  printing  labels  on  tracer  particle 
plots. 

XUM 

(MAP) 

(41) 

■  ~ 

Used,  in  MAP.  Array  has  negative  alphabetic 
characters  foe  maps.  Defined  in  a  DATA 
Statement. 

XX 

B.C. 

( IDXX) 

cm 

Equivalenced  to  X  array  so  as  to  make  X(0) 
available . 

Y 

B.C. 

(JDX) 

cm 

Distance  from  bottom  of  grid  to  top  of  cell 
Equivalenced  to  YY  array  such  that  Y(0) 

-  YY(1). 

YAMC 

B.C. 

(JDX) 

cm-g/sec 

Calculated  and  used  in  PH2.  Axial  momentum 
of  mass  transported  across  left  side  of  cel! 
Equivalenced  to  UL  array.  (See  Appendix  B) 

YCEKT 

*Z(2S) 

cm 

The  tip  of  the  lino- is  described  by: 
y  =  SLOPLN'-x  ♦  YTOPI. 

YMAX 

(SETUP) 

-- 

cm 

=  Y(JMAX). 

YP 

TRACRS 

(KDX4) 

cell 

Y-coordinates  of  passive  cell -centered 
tracer  particles . 

YY 

B.C. 

(JDXX) 

cm 

Equivalenced  to  Y  array  so  as  to  make 

Y(0)  available. 

Z 

B.C. 

Storage  for  most  of  the  input  parameters 
follows  Z,  the  first  word  of  blank  common. 
The  “Z-array"  (the  first  130  words  of  blank 
common)  is  written  on  tape  for  restarting 
problems  (see  Appendix  A). 
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APPENDIX  A 


Z-BLOCK  VARIABLES  LISTED  NUMERICALLY 


See 

Dictionary 

for  meaning 

and  use. 

1. 

PROB 

30. 

NC 

59. 

2. 

CYCLE 

31. 

SLOPLN 

60. 

3. 

DT 

32. 

NRC 

61. 

4. 

NUMSP 

33, 

I  MAX 

62. 

5. 

NFRELP 

34. 

I  MAX  A 

6  3 . 

6. 

N DUMP  7 

35. 

JMAX 

64 . 

7 . 

ICSTOP 

56. 

JMAXA 

65  , 

8. 

PIDY 

37. 

KMAX 

66 . 

9. 

TOP  MU 

38. 

KMAXA 

67. 

10. 

RTMU 

39. 

BOTM 

68 . 

11. 

UN  11 

40. 

BOTMV 

69. 

12. 

NUMREZ 

41. 

NUMSP  T 

70. 

13. 

ET.li 

42. 

MAPS 

71. 

14. 

KUNITR 

43. 

NUMSCA 

72 . 

15. 

IPR 

44. 

PRLIM 

73. 

16. 

PRCNT 

45. 

PRDELT 

74. 

17. 

KUNITW 

46  . 

P REACT 

75. 

18. 

XMAX 

47. 

11 

76. 

19. 

NZ 

48. 

12 

77. 

20. 

NREZ 

49. 

IPCYCL 

78. 

21. 

IGM 

R5  0 . 

TSTOP 

79. 

22. 

NAPRST 

51. 

WFLAGF 

SO  . 

23. 

NAPRND 

52. 

WFLAGL 

81. 

24. 

DM  IN 

53. 

UN  5  3 

82 . 

25. 

YCENT 

54. 

IVARDY 

So  . 

o  1 

26. 

DTNA 

55. 

VT 

S-+ . 

n  7“ 

7  7 

U  I  . 

CVIS 

56. 

N6 

S  b  • 

•  28. 

RTIPI 

57. 

RTM 

S6 . 

29. 

UN  2  9 

5  8. 

RTMV 

S7 . 

UN  5  9 

N10 

Nil 

GAMMA 

TO  PM 

BOTMU 

UN  6  5 

TOP  MV 

NSIDES 

NMAT 

CYCMX 

CYCPH3 

REZFCT 

NT  RAC  R 

NMXCLS 

BBOUND 

UN  7  5 

ECK 

NECYCL 
NTPMX 
UN  7  9 

uvmax 

NTCC 
UN  S  2 
IYARDX 
T 

EMIN 
PM  IN 
INTER 
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88. 

I 

123. 

IEXTX 

89. 

J 

124. 

JEXTY 

90. 

K 

125. 

UN  1 2  5 

91. 

M 

126. 

UN  1 2  6 

92. 

N 

127. 

SSI 

93. 

UN  9  3 

128. 

SS2 

94.  . 

YTOPI 

129. 

UMIN 

95. 

REZ 

130. 

SS4 

96. 

NO  DUMP 

131. 

PRTIME 

97. 

NTCRSV 

132. 

EOR 

CO 

TJMOB 

133. 

EOT 

99. 

XCONI 

134. 

EOB 

100. 

EVAPM 

135. 

EMOR 

101. 

EVAPEN 

136. 

DXF 

102. 

EVAPMU 

137. 

DYF 

103. 

EVAPMV 

138. 

UN  1 3  8 

104. 

EZP1I2 

139. 

-STAB 

105. 

UNI  05 

140. 

XIENRG 

106. 

MATXPL 

141. 

XKENRG 

107. 

XPLRAD 

142. 

XTENRG 

108. 

XPLCNT  . 

143. 

DTFIX 

109. 

XPLVEL 

144. 

DTMIN 

110. 

ROEPS 

145. 

UNI  4  5 

111. 

UN  111 

146. 

EMOT 

112. 

UN  1 1 2 

147. 

JCENTR 

113. 

FINAL 

148. 

RADIUS 

114. 

UN  1 1 4 

149. 

BBAR 

115. 

MBBB 

150. 

EMOB  = 

116. 

MBB 

117. 

118. 

NKDX4 

NADD 

151. 

PK(1) 

119. 

MINX 

152. 

PK(2)  , 

120. 

MAXX 

153. 

PK(3)  : 

121. 

MI  NY 

122. 

MAXY 

0  [Last  card  of  SETUP  input 
deck.  Not  included  when 
restarting  from  tape.] 

should  be  the  same  as  PROB, 
problem  number. 

the  restart  cycle. 

when  =  -1.,  program  will  re¬ 
start  from  tape  and  do  a  "long'1 
EDIT  prints  of  the  pickup  cycle. 
When  =  *2.,  program  will  restart 
from  tape  and  do  a  "short"  EDIT 
print,  of  tlie  pickup  cycle. 
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APPENDIX  B 


VARIABLES  USED  IN  TRANSPORT  ACROSS  CELL  BOUNDARIES 


AMMY  DELEB 


APPENDIX  C 


FLAGS  AND  CONVENTIONS  GOVERNING  INTERFACE  CELLS  AND 
CELLS  CONTAINING  EXPLOSIVES 


Statements 


Meaning 


R1IO  (1  ,M)  =  -1.0 

NVOID  =  NMAT  +  1 

M  =  MFLAG (K) - 100 
RIIO  (NVOID, M)  »  1.0 

MFLAG (K)  =  0 
0  <  MFLAG (K)  <100 
MFLAG (K)  =  2 

MFLAG (K)  >  100 

MFLAG (K)  <  0 


The  M  location  in  mixed  cell  arrays 
(XMASS,  RHO,  SIE,  FRACTP ,  FRACRT, 
VOLMAT,  VOLRT,  VOLTP,  SAMMP,  SAMPY) 
is  not  in  use . 

NVOID,  the  number  of  the  "void" 
package ,  is  one  more  than  the  number 
of  material  packages. 

Cell  K  contains  the  free  surface 
interface. 

Cell  K  is  empty  and  is  not  cut  by 
an  interface,  i.e.,  it  is  a  void. cell. 

Cell  K  is  nonempty  and  pure,  ie .  , 
does  not  contain  an  interface. 

Cell  K  is  completely  inside  package 
2  boundary  and  contains  only  material 
of  package  2. 

Cell  K  is  mixed,  i.e.,  contains  at 
least  one  interface.  M  =  FLAG(K)-100 
gives  location  of  quantities  in  mixed 
cell  arrays  for  cell  K. 

(a)  During  INFACE  or  PI  1 2 ,  cell  K  was 
mixed,  but  is  no  longer  cut  by 

an  interface  and  will  be  reflagged 
at  t lie  end  of  P1I2. 

(b)  During  COT  or  Pill,  cell  K  has  an 
artificial  driving  force  acting 

on  it.  Pressures  are  set  by  APLYPR 
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Statements 


Meaning 


MAT ( 1 )  =  3 


M  =  MFLAG (K) -100 
RIIO  (N,M)  =  0. 


M  =  MFLAG (K) - 100 
RHO  (N ,M)  >  0 
XMASS (N  ,M)  =  0. 


M  =  MFLAG (K) -100 
RHO  (N ,  M)  *  0. 
XMASS (N,M)  >  0. 


MAT (2)  =  20 


NAPRST  =  0 

NAPRST  >  0 
NMAT  =  1 


M ATX PL  >  0 
MAT  (2)  =  20 


The  material  code  number  of  package 
is  3.  The  list  in  EQST  indicates 
that  code  number  3  corresponds  to 
iron. 

The  interface  of  package  N  does  not 
cut  cell  K.  Material  of  package  N 
will  not  be  transported  into  cell  K. 

Package  N  interface . cuts  cell  K, 
but  cell  K  does  not  yet  contain  any 
material  of  package  N. 

Package  N  interface  has  left  cell 
K  and  the  material  of  package  N 
[ XMASS (N,M)]  should  be  completely 
evacuated  on  this  cycle. 

The  material  of  package  .2  is  an 
Ideal  Gas  (MATXPL=0)  or  an  explosive 
(MATXPLrO)  and  is  given  special 
treatment  in  the  pressure  calculation, 
sound  speed  calculation  and  the 
strength  phase  of  the  code. 

No  pressure  boundary  condition  is 
being  used. 

A  pressure  boundary  condition  is 
being  applied  along  a  section  of  the 
free  surface  defined  by  free  surface 
tracers.  [TX(XVOIl),  NAPRST), 

TY  (NVOI 1),  NAPRST)  through 

TX  (N VO  I D ,  NAP RND)  ,  TY  (NVOI  D NAP RXD)  ]  . 

Material  package  2  is  an  explosive 
detonating  with  a  plane  detonation 
front. 
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Statements 


Meaning 


MATXPL  <  0 
MAT (  2)  =  20 


Material  package  2  is  an  explosive 
detonating  with  a  curved  detonation 
front. 
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